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ABSTRACT 

We develop a method to calculate the column density distribution of the 
Lya forest for column densities in the range 10^^'^ — 10^^'^ cm~^. The Zel'dovich 
approximation, with appropriate smoothing, is used to compute the density 
and peculiar velocity fields. The effect of the latter on absorption profiles is 
discussed and it is shown to have little effect on the column density distribution. 
An approximation is introduced in which the column density distribution is 
related to a statistic of density peaks (involving its height and first and second 
derivatives along the line of sight) in real space. We show that the slope of the 
column density distribution is determined by the temperature-density relation as 
well as the power spectrum on scales 2 /iMpc~^ ^ A; ^ 20 hMpc~^. An expression 
relating the three is given. We find very good agreement between the column 
density distribution obtained by applying the Voigt-profile-fitting technique 
to the output of a full hydrodynamic simulation and that obtained using our 
approximate method for a test model. This formalism then is applied to study 
a group of CDM as well as CHDM models. We show that the amplitude of 
the column density distribution depends on the combination of parameters 
{rihh'^yTQ^''^ J^l , which is not well-constrained by independent observations. 
The slope of the distribution, on the other hand, can be used to distinguish 
between different models: those with a smaller amplitude and a steeper slope 
of the power spectrum on small scales give rise to steeper distributions, for the 
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range of column densities we study. Comparison with high resolution Keck data 
is made. 



Subject headings: cosmology: theory — intergalactic medium — quasars: 
absorption lines 



1. Introduction 

There is a long history of theoretical efforts to place the study of the Lya forest within 
the framework of cosmological structure formation theories (Doroshkevich & Shandarin 
1977; Rees 1986; Bond, Szalay & Silk 1988; McGill 1990; Bi, B5rner & Chu 1992). Recent 
work making use of numerical simulations has greatly advanced our understanding in 
this direction ( |Cen et al. 1994] ; Zhang, Anninos & Norman 1995; [Hernquist et al. 1995 ; 



Petitjean, Miicket fc Kates 1995| ; [Miralda-Escude et al. 19961) . (See also Bi, Ge & Fang 



1995. for a linear theory calculation). The emerging picture is that it is possible to account 
for all the observed properties of the Lya forest (with column densities less than about 
10^^ cm~^) by assuming it originates from the small scale structure, including the network 
of filaments, pancakes and mild density fluctuations, which arises naturally in hierarchical 
clustering cosmological models ( [Weinberg et al. 1996|) . 



A commonly used statistic to characterize the forest is its column density distribution, 
the number of absorption lines per unit neutral hydrogen column density per unit redshift 
as a function of column density. Other useful statistics include line-line correlations and the 
distributions of 6- values and equivalent widths ([Murdoch et al. 1986 ; Carswell et al. 1991 



Press, Rybicki fc Schneider 1993[ ; [Crisitani et al. 1995[) . There have also been proposals of 



new statistical tools ([Meiksin and Bouchet 19"95| ; [Miralda-Escude et al. 19"96| ; Pando and[ 



Fang 1996 ). (See Tytler 1992| for a general overview of the statistical issues concerning 



quasar absorption systems.) We focus our attention on the column density distribution in 
the present work. 

One of the most striking features of the observed column density distribution of quasar 
absorption systems is that it can be approximated by a single power law that extends over 
many orders of magnitude. This was emphasized by Tytler (1987), among others, who found 
that in the range 10^^ < A^hi < 10^^ cm~^, the distribution was reasonably well represented 
by a power law, oc N^f with (3 = 1.51 ± 0.02. However, there exists evidence of at least one 
break. It has been demonstrated that there is a deficit of absorption systems somewhere 
in the column density range 10^^ to 10^''cm~^ compared to a power-law extrapolation of 
the distribution from lower column densities ( Parswell et al. 1987] ; [Petitjean et al. 1993[ ; 
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H.U et al. 1995 ; [Giallongo et al. 1996|) . For reasons that have to do with the nature 



of the approximations that we make (Sec. |5.2| ), we focus our attention on absorption 
systems with column densities in the range 10^^'^ < A^hi < lO^^'^cm"^. Hu et al. (1995) 
obtained (3 = 1.46 with a 95% confidence range of (1.37, 1.51) in the column density range 
-|^q12.3 ^ jY^^ ^ 10^^'^ cm~^ . Lu et al. (1996) found the same best-fit (3 for the same range 
of column densities. 

An obvious ultimate goal of recent theoretical work on the Lya forest is to constrain 
theories of structure formation. The natural question is: what determines the normalization 
and slope of the column density distribution? What are the major determining factors, 
in addition to the usual parameters specified by a given cosmological model? To answer 
these questions, another question has to be addressed: what are the analytical and/or 
computational tools necessary to make accurate predictions for the column density 
distribution, given all the required parameters? 

Accordingly, the present work can be divided into three parts, where the tools are 
developed, the factors that influence the column density distribution are analyzed and one 
application to a class of cosmological models is discussed. 



Numerical hydrodynamic simulations (|Cen et al. 1994| ; [Zhang et al. 1995| ; [Hernquist et 



|al. 1995| ; [Miralda-Escude et al. 19961) provide the most obvious tools to study the Lya forest. 



Computational costs, however, prevent one from testing extensively several cosmological 
models. We show in this paper that the Zel'dovich approximation ( [Zerdovich 19701 ), with 



appropriate smoothing, is an efficient and accurate alternative (see also Doroshkevich & 
Shandarin 1977). Our basic assumption is that the part of the Lya forest with column 
densities less than about 10^^'^ cm~^ arises mostly from regions which are slightly overdense 
(overdensity ^ 5) or even underdense and which have not undergone orbit-crossing. The 
Zel'dovich approximation can then be coupled with the equations governing the thermal 
and ionization states of the gas to yield accurate predictions for the density of neutral 
hydrogen and the peculiar velocity as a function of position. Absorption spectra are then 
generated and analyzed. Basic expressions for the absorption optical depth are presented in 
Sec. and the approximations that go into its computation are discussed in Sec. ^ 

Given an absorption spectrum, the column density distribution depends on the method 
of identifying lines and assigning column densities. This is discussed in Sec. [4.1|. We 



investigate the effects of peculiar velocities on the column density distribution, using a 
method described by Miralda-Escude (1996). We find that although peculiar velocities 
can strongly infiuence the shapes of absorption profiles, they play a relatively minor role 
in determining the column density distribution. The various interesting effects of peculiar 
velocities are discussed in Sec. [4.2|. Motivated by this finding, a very different way of 
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assigning column densities is introduced in Sec. in which no absorption spectrum needs 
to be generated. In the absence of pecuhar velocities, there is a one-to-one correspondence 
between density peaks in real space (if they are separated by a distance larger than a 
minimum corresponding to the thermal broadening width) and minima of transmission 
(maxima in absorption) in the observed spectrum. Under such conditions, we can simply 
associate each density peak in real space with an absorption line and assign a column 
density to each based on the height and curvature of the peak. The column density 
distribution is then a statistic of density peaks in real space. We apply this procedure 
(we call it the Density-Peak-Ansatz) to the density field predicted by the truncated 
Zel'dovich approximation and test the result against that of a full hydrodynamic simulation. 
The column density distribution obtained in this way is compared to that obtained 
from the hydrodynamic simulation using the Voigt-profile-fitting-technique, which is the 
line-identification method most widely used. The level of agreement is found to be excellent. 
In Sec. [5.2|, we discuss the range of parameters in which our computed column density 



distribution is expected to be reliable. 

Armed with the right tools, we turn our attention to the second question: what factors 
determine the column density distribution? They can be divided into two categories. One 
has to do with properties of the intergalactic medium, including its temperature, the 
equation of state (or temperature-density relation, which we will use interchangeably; see 
Hui & Gnedin 1996), the ionizing radiation intensity and the baryon density. Uncertainties 
in all of them have to be taken into account before one can meaningfully confront theories 
with observations. We distinguish between the factors that mostly affect the normalization 
of the column density distribution and those that mostly affect its slope. It is found that 
the temperature-density relation (weakly) affects the slope while the rest of the above 
factors influences the normalization. It is also emphasized that the temperature and the 
equation of state depend on the reionization history of the universe (a fuller discussion of 
this point and related topics will be given in a separate paper). The second set of factors 
affecting the column density distribution has to do with the specific cosmological model, 
namely the normalization and shape of the corresponding power spectrum. We study a 
few variants of the Cold Dark Matter (CDM) model in Sec.^ for this purpose. It is found 
that the amplitude and slope of the linear power spectrum on comoving scales of around 
2h Mpc~^ to 2{]h Mpc~^ are the most important factors in determining the slope of the 
column density distribution (the equation of state also has a weak effect on it). Decreasing 
the amplitude and/or steepening the slope of the power spectrum tends to steepen the 
distribution in the column density range about 10^^'^ to lO^'^'^cm"^. We introduce an 
expression relating the slope of the column density distribution to the equation of state and 
properties of the power spectrum on small scales. 
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We then study a class of Cold plus Hot Dark Matter (CHDM) models in Sec. ^ 
making use of the insights gained in Sec. |^ and Sec. |^ The Q^, = 0.2 (density parameter 
in neutrino) models have steeper column density distributions compared to those with 
fijy = 0.1 because they have less power on the relevant scales. In particular, the low Hubble 
constant [Hq = 50 kms~^ Mpc~^) = 0.2 models predict slopes that are steeper than the 
observed one for most of the parameter-space specifying the properties of the intergalactic 
medium. Only for equations of state that are close to isothermal can the two be made 
consistent with each other. We emphasize however that a more detailed comparison between 
the models and observations, taking fully into account instrumental noise and biases of 
the line-identification method(s), is necessary before one can firmly reject any model. We 
conclude in Sec. ^ . 

It is appropriate that we mention here two recent pieces of work along similar lines 
as described above, but using a different dynamical approximation, namely the lognormal 
approximation: Bi & Davidsen (1997) and Gnedin & Hui (1996). The former, in particular, 
contains a very comprehensive and careful analysis of the many different observational 
properties of the Lya forest. One strong point of their analysis is that they tested their 
method using VPFIT, a spectral analysis routine that is commonly used by observers. We 
will discuss the predictions for the column density distribution by the lognormal and the 
Zel'dovich approximations in Sec. ^. 

In our notation, bold faced letters are reserved for three-dimensional vectors. The 
symbols Vpec and x denote the three-dimensional peculiar velocity and comoving position 
while fpec and x are their counterparts along the line of sight of interest. Standard symbols 
are used for cosmological parameters: H for the Hubble constant as a function of z, Hq 
for the Hubble constant today, h for Hq/100 kms^^Mpc^^, Qq for the density parameter 
today, with the subscript b to denote its baryon portion and u its neutrino content. We 
use the symbol h (as distinct from h) to denote the Planck constant in a few places where 
it arises. The term multiple-streaming is reserved for the situation where a single observed 
redshift corresponds to more than one position in real space. We distinguish it from the 
term orbit-crossing, which is commonly used interchangeably with multiple-streaming in 
other contexts. Orbit-crossing refers to the case where a single position has more than one 
velocity. 
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2. Cosmological Lyman-Alpha Absorption in a Fluctuating Medium: Basic 

Results 

A photon emitted with energy higher than 10.196 eV (wavelength of 1216 A) by a 
distant quasar is continuously redshifted as it travels through the intergalactic medium until 
it reaches the observer. At some intermediate point, the photon is redshifted to around 
1216 A in the rest frame of the intervening medium, which may contain neutral hydrogen. 
It can then excite the Lya transition and be absorbed. Let us consider a particular line of 
sight from the observer to the quasar. The optical depth r (the probability of transmission 
is given by e~^) of a photon at a given observed frequency is given by: 

t{uo) = / nmaa—— , (1) 

JxA i- + Z 

where x is the comoving radial coordinate of some intermediate point along the line of sight, 
z is the redshift and nni is the proper number density of neutral hydrogen at that point. 
The limits of integration, xa and xb, are the comoving positions of the observer and the 
quasar. The Lja absorption cross-section is denoted by cTq,. It is a function of the frequency 
of the photon with respect to the rest frame of the intervening hydrogen at position x. Let 
us call this frequency u. The cross-section is peaked when u is equal to the Lya frequency. 

The frequency u is related to the observed frequency Uo by: 

+ , (2) 

where Vp^c is the peculiar velocity along the line of sight at position x and 1 -|- 2; is the 
redshift factor due to the uniform Hubble expansion alone at the same position. The 
peculiar velocity of the observer, which merely displaces the whole spectrum by a constant 
amount (independent of x), is ignored. The quantity Vpec/c, where c is the speed of light, is 
much smaller than 1. 

It proves convenient for later discussion to expand z around some mean redshift of 
interest z, which could be the redshift of a simulation output or the average redshift of an 
observed spectrum with limited redshift range. Using dx — cdt/a, where a is the Hubble 
scale factor and t is the proper time, it can be shown that 

i/ = i/o(l + ^) fl + , U = -^^{x - X) + Vpec{x) , (3) 

V c/ 1 + z 

where x is the position at which the redshift due to Hubble expansion coincides exactly 
with z. The Hubble constant at the same redshift is denoted by H. We assume the range 
of X is small enough so that u/c<^ 1. The convention a — 1 today is adopted. 
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The velocity coordinate u defined above contains contributions from botli tlie Hubble 
expansion and the peculiar motion. Without peculiar motion, u increases monotonically 
with X and is in fact linear in x. Peculiar velocities destroy the linear relation and could 
give rise to situations where a given u corresponds to more than one position x. It implies 
that a photon of a given observed frequency z/q can have the same rest-frame frequency v at 
more than one place in its trajectory from the quasar to the observer. We reserve the term 
multiple-streaming to this situation and distinguish it from orbit-crossing where a given 
X carries more than one Vp^c or u. We will return to the subject of multiple-streaming in 



Sec. and that of orbit-crossing in Sec. 3.1 



We define one more velocity coordinate Mq, which is related to the observed frequency 
i/o by: 

-0 = ^(1-^) (4) 

1 + z V c / 

where Va is the Lya frequency. The velocity coordinate Uo is simply equal to u when u 
coincides exactly with z/q, (this can be seen by comparing eq. and bearing in mind 
that u/c and Uo/c are both assumed to be much less than 1). 

With the definitions in place, we change the variable from a; to u in equation (|I|), which 
results in the following expression for r, now a function of Uo'- 



T[Uc_ 



E 



Ua 1 + ^ 



du ' ^ 



dx 



aadu , a„ = a^or^e-^^-"")'/"' . (5) 

by/TT 



The summation refers to a sum over multiple-streams (all the x's within the range xa — xb 
that corresponds to a given u), and nni, z and \du/dx\~^ are now functions of u. The limits 
of integration ua and ub are the velocity coordinates corresponding to the positions xa and 
Xb (assuming no orbit-crossing so that each x carries one u). Note that in practice, only 
a limited range of u contributes to r for a limited range of Uo so that one can replace the 
redshift z with z. The same is also true for equation (|l|). 

The Lya cross-section is expressed as a function of u — Uo- The constant aao is equal to 
the combination of fundamental physical constants 0.4167rg^/(meCz/a), where q is the charge 
of an electron and mg is its mass. It is about 4.5 x 10~^^cm^. The parameter b is equal to 



2kBT/mp where fc^ is the Boltzmann constant, nip is the mass of a proton and T is the 
temperature of the gas at the velocity coordinate u. 

The form of the line profile function above takes into account thermal broadening but 
ignores the natural line width. A more general profile function involves a convolution of the 
two, resulting in the Voigt profile ( [Rybicki fc Lightman 19791) . However, the Voigt profiles 



are accurately thermal profiles for column density less than about 10 cm . The reader is 
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referred to Spitzer (1978) and Press and Rybicki (1993) for discussions of curve of growth 
analysis. 

Note also that it is sometimes assumed b contains a component due to turbulent 
motion. We do not include it explicitly in our formalism. Bulk motion, on the other hand, 
is accounted for by fpec or u. 

Let us consider two different limits of equation (^). 

Suppose there is a high local maximum in nm\du/dx\~^ at some u = Umax with width 
in velocity space much smaller than the thermal width b. Then one can take the line profile 
function associated with cTq, out of the integral in equation (^ because niii\du / dx\~^ varies 
much more rapidly than the thermal profile: 

r(wo) = ( I nHi(x)^)a,o^e-(-— , (6) 



max 



where the variable of integration has been changed back from u to x. The equation holds 
when Uo is close enough to Umax- The integral is over the local maximum, assuming that the 
amount of neutral hydrogen away from the maximum does not cause significant absorption 
(until another maximum is encountered). One then sees an absorption line with a Gaussian 
profile in optical depth. While the width of the line tells us about 6, which is proportional 
to the square root of the temperature, the depth of the line provides information about 
both b and the column density, which is the integral inside the first pair of brackets on the 
right hand side. Let us call this the narrow- maximum-limit. 

Consider another limit of the integral (eq. [^) in which n^uldu / dx\~^ varies slowly with 
u. Suppose the scale of variation is much larger than the thermal width. In this case, one 
can leave the line profile function inside the integral but take the rest outside:[] 



( \ ""-HI 



+ z 



du ' ^ 



dx 



caao ■ (7) 



The velocity dependent terms on the right hand side are evaluated at Uq. The profile 
function has been integrated out. 

In the above limit, r does not necessarily have the thermal profile around its maxima. 
We will call this the broad- maximum-limit. An extreme example is that of a homogeneous 
medium, which gives rise to featureless and uniform absorption ( |Gunn fc Peterson 19651 ). 



''The expression is not valid at velocity caustics, where du/dx vanishes. Further discussion on velocity 



caustics can be found in Sec. 4.2 
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Conventional analysis of quasar spectra involves identifying those parts of the spectra 
that are due to the Lya absorption and fitting them with superpositions of the Voigt 
profiles (of which the thermal profiles are a subset) until the residual signal is consistent 
with noise. This technique was motivated by the picture of the intergalactic medium as 
consisting of a smooth component which causes relatively little absorption and a set of 
clouds that satisfy the narrow-maximum-limit. For each cloud, the best-fit Voigt profile 
then yields its temperature and column density according to equation (H). 

However, it is clear that for a general fluctuating medium, not all maxima in r 
necessarily satisfy the conditions leading to equation (^. In fact, according to most 
structure formation theories, there invariably exist fluctuations in the intergalactic medium 
on scales larger than the thermal width. In the broad- maximum-limit, the shape of a local 
maximum in optical depth is determined by the distributions of whi and \du/dx\ around it. 
Each maximum in r can still be identified as an absorption line and one can even apply 
standard techniques and try to fit its shape with superposition of the Voigt profiles. Given 
the best-fit Voigt profiles, one can assign a 6- value (width of the profile) and a column density 
to each profile but it is no longer true, for instance, that the 6-value thus obtained is equal 
to \j2kBT /nip (eq. j^). A reasonable question to ask is whether there are other practical 
methods of identifying absorption lines and assigning column densities without assuming 
every absorption line consists of a superposition of the Voigt profiles. This will be discussed 
in Sec. |4.1| . It should be borne in mind, however, that all existing observational data on the 
column density distribution are obtained using the Voigt-profile-fitting-technique, so for the 
purpose of comparing theory with these observations, it is important the line-identification 
algorithm one uses gives results consistent with the profile- fitting-technique. 

In general, there are regions of high density and limited extent, galaxies for instance, 
which give rise to absorption profiles well approximated by the narrow- maximum-limit, 
but there are also regions in the intergalactic medium with gentle fluctuations where the 
broad-maximum-limit holds. Then there are those places where neither limit applies, in 
which full integration of equation (||) has to be carried out to compute the optical 

depth. To do so, one needs to know how the neutral hydrogen density, peculiar velocity and 
temperature vary with spatial position. This is the subject of the next section. In any case, 
the above discussion should make it clear that the quasar absorption spectrum contains a 
wealth of information on the intergalactic medium. 
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3. Ingredients for Generating Quasar Absorption Spectra 



There are four quantities that go into the computation of the optical depth: 
temperature, peculiar velocity, overdensity and neutral fraction. That the temperature and 
peculiar velocity are important should be obvious from the expression for the absorption 
cross-section in equation @. The temperature determines the extent of thermal broadening 
(6) and the peculiar velocity changes the frequency of the photon in the fluid rest-frame 
(eq. [^). Let us define carefully what we mean by the other two quantities, the overdensity 
and the neutral fraction. 

Suppose nH(x) is the total proper number density of all hydrogen species at position x 
and nn is its spatial average. The overdensity 6^, which describes the variation in space of 
rzH(x), satisfies: 



In the first expression, (5b is defined as the number overdensity of hydrogen. In the 
second expression , we equate 6^ with the mass overdensity of baryons {pb is the proper 
mass density of baryons and pb is its mean), which is an excellent approximation for our 
application because there is no significant conversion of hydrogen into other elements, nor 
is there any interaction that could cause the spatial distribution of hydrogen to deviate 
significantly from that of other types of baryons. 

What the Lya absorption directly probes is not the total hydrogen density but its 
neutral component. The neutral fraction Xhi is defined by the following relation: 



where tt-hi is the proper number density of neutral hydrogen as a function of space. The 
neutral fraction is determined by the balance between recombination and ionization, the 
rates of which are dictated by the temperature and radiation intensity respectively. 

In general, all four quantities, overdensity ^b, peculiar velocity fpec, temperature T and 
neutral fraction Xhi, are functions of position. In the next two sub-sections, we discuss 
first how to determine the spatial distributions of 5b and Upec, and second how T and Xhi 
vary with position through their dependence on 5b. AH quantities are evaluated at z = 3. 
Although most of the material in this section is standard textbook fare, it consists of a 
somewhat unusual combination of methods, so it is worth going through the basic equations 
and stating our assumptions carefully. 






(9) 
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3.1. The Zel'dovich Approximation 



In cosmological models where dark matter (a term we use interchangeably with 
non-interacting matter) dominates the mass density of the universe, 6^ as defined in 
equation (||) coincides with the dark matter overdensity 5dm on large scales. We define (5dm 
in an analogous manner as before (eq. |^): 



where pdm is the mass density of dark matter at position x and pdm is its mean. The 
equality 5b = 5dm is equivalent to the statement that the hydrogen density (which we 
assume is simply proportional to the baryon density) varies with position in the same 
manner as the dark matter density does. This is true on large scales where gas pressure 
is insignificant compared to the gravitational attraction of the dark matter, provided 
the baryons and dark matter start out having the same spatial distribution, which is 
approximately true for most popular cosmological models. Moreover, with no significant 
interaction that distinguishes between the two on large scales, the baryons and dark matter 
share the same peculiar velocity field. On small scales, however, gas pressure can cause 
the spatial distributions of baryons and dark matter and their velocities to differ. We will 
return to this point below. 

Hence, on sufficiently large scales (how large is large, an obviously important question, 
will be addressed later), it is adequate to know the overdensity and peculiar velocity of the 
dark matter. The Zel'dovich approximation ( [Zerdovich 19701 ) is a well-tested approximation 
to compute the density and velocity distributions of dark matter in the mildly nonlinear 
regime (5dm ~ 5) before orbit-crossing. The reader is referred to the article by Shandarin 
and Zel'dovich (1989) for a comprehensive review (see also Hui and Bertschinger 1996 for 
an alternative formulation of the approximation). 

The starting point of the Zel'dovich approximation is the following equation for the 
displacement of a given mass element or particle: 



The coordinate q is the initial position of the mass element and x is its comoving position 
as a function of time. The displacement is then D_(_(t)VqV'(q)- Its time dependent part 
Dj^{t) is the linear growth factor (Peebles 1980), which, for a universe with critical matter 
density, can be equated with a, the Hubble scale factor. The time independent function 
VqV'(q) is determined by initial conditions. Growing mode initial conditions dictate that it 
is curl-free, hence its form as the gradient of the potential ijj (Vq is the spatial gradient in q 
space) . 



Pdm(x) = Pdm[1 + '5dm(x)] 



(10) 



x(q,t) = q + Z}+(t)Vq^(q) 



(11) 
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Expressions for the peculiar velocity and overdensity follow immediately from equation 

(0): 

r 

Vpec = aD+Wqip , 1 + 5dm = det"-^ ■ ^ '^'^ 



(12) 



dqidqj_ 

The dot in the first expression denotes differentiation with respect to proper time t. The 
peculiar velocity is defined by Vpec = ad:K/dt. The second expression follows from mass 
conservation i.e. (1 + Sj)M)d^x = d^q. The right hand side of the second expression is simply 
the Jacobian of the q-x transformation matrix. 

The function ^(q) contains all the information on the specific cosmological model one 
chooses to study. For the cosmological models we study in this paper, it is a Gaussian 
random field in q space. Suppose is its Fourier counterpart defined hj ip = J d^kip e^^''^. 
The two-point correlation of ip is related to the commonly used power spectrum P by 

(^(k)V^*(k')) = k-^P{k)6^{k - k') , (13) 

where P is related to the root-mean-squared (rms) linear overdensity fluctuation by 

/■oo 

{5^)=Dl{t) / A7cP{k)edk. (14) 
Jo 

Note that Dl = {1 + z)-^ for a universe with critical matter density, choosing = 1 
today. 

To produce a realization of the density and velocity fields for a given cosmological 
model, we employ the following procedure: first, we use the corresponding power spectrum 
to generate the Gaussian random field ip{q) on a grid; second, we displace particles from 
their initial grid positions (q) according to equation (0) for -D+(t) corresponding to ^ = 3; 
third, a peculiar velocity is assigned to each particle according to the first expression in 
equation (p!2D; finally, we use the TSC (Triangular-Shaped density Cloud) scheme ([Hockney 



fc Eastwood 1988D to interpolate the particle positions and velocities to become momentum 



and mass densities on the grid and divide one by the other to obtain the velocity itself. 
The interpolation to obtain mass density is our way of enforcing mass conservation, as is 
expressed in the second expression of equation (O). In the last procedure, we smooth the 



momentum and mass density fields over a small number of grid cells (in fact, we use one 
and have checked that the precise number is not important as long as it is small) before 
performing the division to obtain the velocity field so that we have well-defined velocities 
even in places with zero density after the TSC interpolation (|Kofman et al. 1994|) . Any line 
of sight can then be chosen through the simulation box and the above set of steps gives the 
overdensity and peculiar velocity (in fact only the component parallel to the line of sight is 
needed) at each grid point on it. 
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The procedure just outlined is very efficient because there is no need to integrate 
any equation of motion. One simply multiplies the displacement of each particle by an 
appropriate factor of D^{t). However, the ffist step of the procedure has to be slightly 
modified to address two problems. 

The ffist one is orbit-crossing. The Zel'dovich approximation is known to predict too 
rapid growth of the thickness of the post-collapse pancake (Shandarin and Zel'dovich 1989). 
A number of cures have been proposed ( [Kofman, Pogosyan, fc Shandarin 1990] ; [Matarrese et| 
al. 1992| ; Prainerd, Scherrer, fc Villumsen 1993| ; [Bagla fc Padmanabhan 1994| ) but the one 



that consistently gives good agreement with N-body simulations is the truncated Zel'dovich 
approximation (Kofman 1991; [Coles, Melott fc Shandarin 1993|; [IVlelott, Buchert fc WeiJj 



The basic idea is to smooth the initial power spectrum on small scales so that the 
amount of orbit-crossing that might have occurred by the time of interest is not significant 
enough to destroy the accuracy of the Zel'dovich approximation on large scales, where the 
fluctuations are still mildly nonlinear. The initial power spectrum P{k) is multiplied by a 
Gaussian smoothing kernel of the form e~'^^/'^s, before it is used to generate the Zel'dovich 
displacement field (eq. [0). (This is equivalent to smoothing the initial density field (5(x) 
through the following convolution: (27rA;^^)~°'^ / 5(x')e~^'^sl^~^'P/2)(^3^/ rjj^g smoothing 
wavenumber A;s is chosen according to the following prescription 

ks = l.SfcNL, where 1 = Dl{t) / A7cP{k)k'^dk . (15) 

Jo 

Note that P{k) above is the initial power spectrum before any smoothing. The 
proportionality constant between kg and /cnl actually depends somewhat on the power 
spectrum, with more smoothing (smaller fcs) required for models that have relatively more 
power on small scales ( [Melott 1994| ). The choice above has been shown to give good 
agreement with N-body simulations for CDM models (Melott et al. 1995). We will see 
that for those CHDM models with relatively little power on small scales, the precise value 
of ks is not important. The procedure described above is commonly called the truncated 
Zel'dovich approximation. 

The second problem is one we have pointed out before, namely that 5dm is not 
necessarily equal to (5b (which is what we are interested in ultimately) on small scales. 
In linear theory, it is possible to show that for a dark matter dominated universe, 
the Fourier components of the two quantities obey 5b (^) = 5dm (fc) when -C /cj and 
Sb{k) = A;j5dm(^)/^^ when k ^ kj. Under some restrictive assumptions (see Appendix A), 
it can be shown that 5hik) = 5dm(^)(1 + k'^/kj)^^- The quantity kj^ is known as the Jeans 
scale and is defined by: 



^^^^ (16) 



\ Aira'^GiipBM 
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where kB is the Boltzmann constant, f is the spatially averaged temperature of the gas, 
/i is the mean mass per particle (for fully ionized gas with primordial abundances, it is 
about O-QrUp where nip is the mass of the proton) and 7 describes the relation between the 
temperature T (the actual, not average, value) and 1 + 6\y: T oc (1 + 6hy~^. Note that 7 
does not necessarily equal 5/3 unless the gas behaves adiabatically. The proofs of the above 
assertions can be found in the Appendix A (see also Bi et al. 1992; Peebles 1993| ). It is 



sufficient to note here that in the linear regime, the baryon density field is smoother than 
that of the dark matter on small scales due to the effect of gas pressure. 

Now, the above relations between 5b and 5dm hold only in the linear regime when both 
quantities are small. To take into account the effect of gas pressure in the mildly nonlinear 
regime, one possibility is to smooth the initial power spectrum by a factor of (1 + k'^/k'j)^'^ 
before generating the displacement field, similar to what is done in the case of the truncated 
Zel'dovich approximation. This method was used by Reisenegger & Miralda-Escude (1995) 
to study the fiuctuating Gunn- Peterson effect. In practice, we smooth the initial power 
spectrum by a Gaussian kernel e~^^^^^^ and find that the two ways of smoothing give very 
similar column density distributions. 

To give an idea of scale, for 7 = 1.5, T = lO'^K and a universe at critical density, fcj 
is equal to 16.8 /iMpc"^. It turns out that for all models considered in this paper except 
the VL^ = 0.2 CHDM models, the truncation scales fcg ^ according to equation (|15]) are 
larger than kj^ (eq. [|l^), for reasonable ranges of temperature and 7. For these models, 
it is unnecessary to smooth the initial power spectrum over the Jeans scale because the 
truncated Zel'dovich approximation already prescribes more smoothing. The opposite 
is true for the Qi, = 0.2 CHDM models. In fact, the amount of small scale power is so 
insignificant for these models that the precise scale of smoothing does not affect the column 
density distribution for column densities of interest (Sec. Orbit-crossing is probably not 
very severe for this class of models. Hence, uncertainty in the Jeans smoothing scale due to 
uncertainties in the temperature and equation of state of the intergalactic medium is not a 
concern. 



To sum up, we smooth the initial power spectrum on the scale of fcg ^ (eq. |T5 



or 



kj^ (eq. Il^]), depending on which is bigger, before it is used to generate the displacement 
field (eq. |]13[) (except for a few test-cases discussed in Sec. 0). The rest of the procedure 
to obtain the overdensity 61, and peculiar velocity fpec on a grid follows as before. The 
implicit assumption underlying the whole procedure is that fiuctuations on scales smaller 
than the smoothing scale do not contribute significantly to the number of absorption lines 
at our column densities of interest, about 10^^'^ — 10^^'^ cm~^. The upper limit is related 
to the maximum overdensity {Sb around 5) beyond which the Zel'dovich approximation 
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is not expected to be reliable and the lower limit is set by our resolution (see Sec. ^ 
for more details). Note that while the Jeans scale smoothing is meant to capture the 
actual smoothing of the density field by gas pressure, the nonlinear scale smoothing is an 
approximation technique to avoid the problem of orbit-crossing. As such, the validity of the 
latter in the present application has to be checked. We show in Sec. |] a comparison between 
the column density distribution computed using the approximate method described here 
and that using a full hydrodynamic simulation. The level of agreement lends support to our 
assumption. Another consistency check is to see if shock-heating is important for regions 
with overdensities (or underdensities) associated with the above range of column densities. 
A plot of density versus temperature like Fig. 2 in Weinberg et al. (1996) shows that 
shock-heating, and by extension orbit-crossing, is not important for regions of underdensity 
or low overdensity. 



3.2. The Thermal and Ionization State 

Given the evolution of 6\, predicted by the Zel'dovich approximation, it is possible to 
integrate evolution equations for T and for Xhi as well as the abundance of other species to 
obtain their relations with 6^,. Q 

Details of the computation will be given in a separate paper (Hui and Gnedin 1996). 
A brief discussion can be found in Appendix B of this paper. We summarize the main 
relevant conclusions here. 

First, ionization equilibrium is maintained at high accuracy except during the period 
of initial reionization. Ionization equilibrium implies that the neutral hydrogen fraction (eq. 
IP]) satisfies 

X,, . 1.6 X 10-« f ^) (^] " (1 + 4) C-±^)\ (17) 

where we have adopted the approximate form of the recombination coefficient of hydrogen: 
4.29 X 10~^^(T/10^K)~'''^cm^s~^, which is sufficient for our purpose (see Hui & Gnedin 1996 
for a more accurate analytical fit). The quantity Jhi is a measure of the radiation intensity 



^ Since these equations are local, in the sense that each mass element evolves independently of the others, 
there is actually no need to generate a full three-dimensional realization for the purpose of studying the 
thermal and ionization evolution. A simpler approach is to generate a set of eigenvalues of the deformation 
matrix d^ip/dqidqj according to the prescription of Doroshkevich (1970) and determine the density evolution 
through the second part of equation (p^. 
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defined as follows (analogous to the definition in Miralda-Escude et al. 1996 but differs by 
a factor of 10~^^ ergs Hz~^ s"^^ cm~^ ster~^): 

Jm = ! (10-^^ ergs Hz"^ s'^ cm'^ ster'^)"^ , (18) 

where Ji, is the specific intensity as a function of frequency u in the units given above, 
h is the Planck constant, hz/Hi is 13.6 eV and am is the ionization cross-section. The 
photoionization rate is simply equal to 4 x 10~^^Jhis~^. 



Observations indicate that Jhi is between about 0.1 and 2.0 for z = 2 — 4 ( |Batjlik. 



Duncan fc Ostriker 198^ ; |Lu et al. 199T| ; [Bechtold 1991 IGiallongo et al. 1996| ; |Cooke 



[Espey fc Carswell 1997] ). A perhaps more common way of characterizing the radiation 
intensity is to quote its value, often referred to as J912, a.t u = um or at wavelength 912A, 
in units of ergsHz"^ s~^ cm~^ ster~^. The relation between J912 and Jhi depends on the 
spectrum. A good approximation for reasonable slopes of the spectrum right above z/hi 
{Ju oc z/~™ for m between 1 and 1.5) is Jhi = 0.7J9i2/10~^^. 

Second, we find that 

T = To(l + 5b)^-\ (19) 

where Tq is not dependent on position, is a good approximation for overdensities of interest, 
(5b ~ 5, with a little flattening at the low end ((5b close to 0.1) for some reionization scenarios. 
We will call this our equation of state. Note that this implies that the spatial dependence 
of T (and by extension Xhi) is completely determined by that of (5^, which is true for 
unshocked gas. Similar power law relations between the overdensity and the temperature 
can be seen in Fig. 2 of Weinberg et al. (1996) for low overdensity. 

Third, both Tq and 7 depend on the reionization history, the reasonable ranges being 
1.2 < 7 < 1.7 and 3000 K < Tq < 30000 K at 2 = 3. It is shown in Hui & Gnedin (1996) 
that 1.3 < 7 < 1.62 at ^ = 3 if the universe reionizes before z = 5, assuming uniform 
radiation field. We allow for a larger range here. 

Combining equations ([T9|) and (p!?]), it can be seen that the neutral hydrogen fraction 
is proportional to (1 + 6, 



,1-0.7(7-1) 



In conclusion to Sec. ^ we have outlined a procedure to use the Zel'dovich 
approximation, with appropriate initial smoothing, to produce a realization of (5b and Vpec as 
a function of position, and we have also shown how the relations between T, Xhi and 6^ can 
be obtained (eq. and [|17l). All of them enter into the calculation of the optical depth 



T (eq. m or [Q). We can compute e called the transmission, which is the ratio of the 
observed to the emitted intensities. Observationally, its measurement requires knowledge of 
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the quasar emission spectrum. Moreover, one must carefully choose the range of frequencies 
to consider if one is to limit the source of absorption to that due to the Lya transition. For 
a discussion of these issues, the reader is referred to Press, Rybicki and Schneider (1993). 
To produce a realistic spectrum, one should also add noise and convolve the transmission 
with a window function to mimic instrumental resolution. This is important for a detailed 
comparison between theories and observations, which we will defer to latter work. Our x 
space grid cells, depending on the particular simulation, have sizes ranging from 0.028 Mpc 
to 0.075 Mpc comoving. Note that the true resolution in velocity space is not uniform 
because peculiar velocity varies from one place to another. Without peculiar velocity, the 
above grid cell sizes correspond to velocity cells of 2.8 — 7.5 kms~^, for h = 0.5 at z = 3 
(eq. 0]). The true velocity resolution is probably a little worse than that. As a comparison, 
high quality Keck Telescope data have a Full- Width-Half-Maximum of about 7kms~^ and 
signal to noise per pixel of the order of 30 or higher ([Hu et al. 1995|; |Lu et al. 1996|). 



4. The Peculiar Velocity: its Effects on Line Shapes and the Column Density 

Distribution 

We show in this section that while the peculiar velocity plays an important role in 
determining the absorption profiles, its effect on the column density distribution is minor. 
The procedures to obtain the column density distribution are discussed first. 



4.1. Line Identification and the Column Density Distribution 



Fig. and Fig. ^ show the velocity, density and transmission (e~^) along two lines of 
sight for a (jg = 0.7 CDM simulation, with h = 0.5 (see Table p. The significance of the 
dashed transmission profile will be explained in the next sub-section. The thermal and 
ionization parameters are described in the caption of Fig. |I|. The truncation scale ks (eq. 
||15|| ) is 2.3 Mpc~^. The transfer function is taken from Ma (1996). We find that using 
instead the transfer function of Bardeen et al. (1986) makes almost no difference to the 
resulting column density distribution, for the range of column densities considered. 



The first thing to note is that for the given parameters. 



b = 



\ 



2kBTo 



13 kms-^(l + (5b)^/^ 



(20) 



This might seem to be too small because the observed lower limit of the 6-value is about 
15 — 20 kms^^ (Pu et al. 1995| ; P!jU et al. 1996|) . There are two points to be made. First, 
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To is the temperature of the gas at Sb = 0. The observed fe-values might originate from 
regions at higher 6h. It is true though that the 1/4 power of 1 + 5f, does not help very much. 
Second, a distinction should be made between the observed 6-value and the b defined above. 
The observed 6-value is obtained by fitting the quasar spectrum with superpositions of the 
Voigt profiles. Each Voigt profile yields a column density and a 6-value. All the density 
peaks that give rise to absorption troughs in Fig. ^ and Fig. ^ have velocity widths larger 
than or comparable to the small thermal width defined in equation (|20|) . Therefore the 
narrow- maximum-limit (eq. 0) does not apply and the absorption troughs do not exactly 
have the Voigt profile shapes. The 6-value obtained from the best-fit Voigt profile of a 
given absorption trough does not necessarily correspond to the thermal width in equation 
(pop. It should also be emphasized that the recent hydrodynamic simulations of the Lya 
forest, which have been so successful in accounting for a lot of its observed properties, have 
similarly low temperatures (see for instance Weinberg et al. 1996 0). 

One might wonder if there exists an alternative spectrum reduction method where 
the Voigt profile is not assumed to be the fundamental shape of absorption troughs, and 
for such a method, how the column density is assigned to each trough. The Voigt-profile- 
fitting-technique is nonetheless very important because it is how all existing observational 
data on the column density distribution are obtained. 

An alternative line identification algorithm was proposed by Miralda-Escude et al. 
(1996) and was also used by Hernquist et al. (1995). First a transmission (e~^) threshold is 
chosen. Any part of the spectrum that is continuously below the threshold is identified as 
an absorption line. The column density associated with it is defined by 

Nni^ — [, riu.)—, (21) 

0" Q-o J line C 

where aao is defined after equation (||). The limits of integration are taken to be over the 
absorption line, i.e. where the transmission is continuously below the threshold. Note that 
if the narrow-maximum-limit or the thin cloud assumption were to hold, equation (|]) can 
be substituted into equation (|2T|) to show that Nm does correspond to / riuidx/^l + z), 
assuming the threshold is high enough so that most of the Voigt profile is included in the 
definition of the absorption line. 

Let us call the above procedure the Threshold- Algorithm. We show in Fig. |^ the 



^Their output is at redshift of 2 and so naturally they have a lower temperature. In general, the 
temperature Tq is dependent on the reionization history of the universe: crudely speaking, the closer the 
epoch of reionization is to the epoch of observation (z — 3 in our case), the higher the temperature. Assuming 
reionization occurs before a redshift of 5, say, puts an upper bound on To (Hui and Gncdin 1996). 
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column density distribution computed according to the algorithm (Crosses). The symbol 
(PNi^ya/ dN-iii/ dz denotes the number of absorption lines per unit column density per unit 
redshift. The reason for the chosen range of column densities will be given in Sec. |^. The 
transmission is chosen to be at the mean transmission of 0.89. 

The Threshold-Algorithm has the tendency to underestimate the number of absorption 
lines compared to the Voigt-profile-fitting-technique. One reason is that it does not deblend. 
In other words, a given absorption line according to the Threshold-Algorithm may contain 
more than one minimum in transmission. Such an absorption line would be broken up to a 
few lines if the Voigt-profile-fitting-technique is employed. To demonstrate this effect, we 
modify the Threshold-Algorithm: for each (parent) absorption line identified, we break it up 
into individual components (children) where each component is bordered by local maxima 
in the transmission within the confines of the parent. The column density for each child 
component is defined similarly as in equation (^Tf ) and the limits of integration are taken to 
be the boundaries of each component. We will call it the Threshold-Deblending- Algorithm. 

The resulting column density distribution is denoted by the square symbols in Fig. ^ 
for the transmission threshold of 0.89. One can see that indeed the number of lines of low 
column densities go up. We should emphasize however that the Threshold-Deblending- 
Algorithm cannot be used to analyze observational data without modifications because in 
real life, noise creates local transmission maxima within any parent absorption line.|^ 

For now, the Threshold-Algorithm is adopted as a simple way to identify lines and 
assign column densities, which we will use to study the effects of the peculiar velocity on 
the column density distribution. 

4.2. The Role of Peculiar Velocities 

The following experiment is performed to investigate the importance of peculiar 
velocities. We generate absorption spectra and compute the column density distribution 
using the same density field as that used to produce the solid curves in Fig. |I], Fig. |^ and 
the crosses in Fig. |^ but we set all peculiar velocities to zero. 

Let us first examine some examples of the absorption spectra. The dashed curves in 
Fig. |l] and Fig. || are the resulting spectra after putting all peculiar velocities to zero. 



^"^In fact, numerical noise can also have the same effect. We check it by defining local maxima in two ways: 
local maxima over three cells and local maxima over five cells with the slope on either side of the maxima 
not changing signs, ft turns out the resulting column density distributions are almost the same. 
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A comparison of the dashed absorption spectrum with its sohd counterpart in each 
figure shows that the pecuhar velocities play an important role in determining the shapes 
of absorption lines. Without peculiar velocities, the shapes of absorption troughs mirror 
closely (perhaps a little thermally-broadened compared to) those of the density peaks while 
with nonzero peculiar velocities, the absorption troughs can have quite different shapes 
from the underlying density field. Peculiar velocities can add or erase structures. An 
example of the former can be found in the pair of density peaks around x = 9 Mpc and their 
corresponding absorption profiles in Fig. 0. An example of the latter can be found in the 
density peak(s) around x = 7 Mpc and the corresponding absorption trough(s) in Fig. |[ 

Broadly speaking, the effects of peculiar velocities on absorption spectra fall into three 
categories. They are distinguished by the value of du/dx {u is defined in eq. |0]). First, 
there are regions in space where the peculiar velocity gradient is small so that du/dx is 
almost equal to its Hubble value H / (1 + z) (eq. [^). An example is the density peak around 
X = 2.2 Mpc in Fig. |^. The peculiar velocity shifts the position of the associated absorption 
trough but does not affect its shape. 

Second, there are places where the peculiar velocity gradient is opposite in sign and 
comparable in magnitude to H/ (1 + z), in which case \du/dx\ becomes very small. Suppose 
also that \d'^u/dx'^\ is small. The implication is that a small range in u corresponds to a 
relatively large range in x. See for instance the density peak(s) around x = 7 Mpc in Fig. H, 
which is a really broad structure in x space but is relatively narrow in u space if peculiar 
velocities are not put to zero. The small \du/dx\ or the converging peculiar velocity flow 
around it helps to produce a narrow absorption trough (second panel from the top in 
Fig. Contrast it with the corresponding absorption feature in the top panel of the same 
figure, where peculiar velocities are set to zero. The limiting case where \du/dx\ exactly 
vanishes is called a velocity caustic ( [McGill 19901 ). 

Third, there are regions where the peculiar velocity gradient dominates in such a way 
that du/dx is negative and \du/dx\ is not small. An example can be found around the pair 
of density peaks at x = 9 Mpc in Fig. ||. This is where multiple-streaming occurs. A given 
range in u corresponds to disjoint pieces in x space. As a result, the shapes of the associated 
absorption troughs are significantly different from those of the underlying density peaks. 

The three categories can be shown to correspond to the three evolutionary stages of a 
pancake collapsing along the line of sight (McGill 1990). Restricting equations (|lID and 
(|l^) to one dimension and putting it into equation (^, one obtains: 
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where we have assumed a universe of critical matter density so that = a. Restricting 
equation ( ]T2| ) to one dimension, it can also be shown that 

1 + 5b = (^1 + , (23) 

where we have replaced 5dm by 6^ assuming the appropriate initial smoothing has been 
carried out as indicated in Sec. As a grows, it can be seen that du/dx goes through the 
three different regimes outlined above for negative d^ip/dq^. At the velocity caustic where 
du/dx = 0, it can be shown that 6^ = 1 QMcGill 1990|) . This conclusion does not hold in 
general of course because pancakes can collapse in directions different from the line of sight. 
But it is true that velocity caustics are often found in regions of slight overdensities. 

In principle, at a velocity caustic, an absorption line can arise even without any 
variation in the density field at all (eq. p). In practice, one expects that converging 
peculiar velocity flows are accompanied by density peaks. This is consistent with the few 
examples we have seen. 

Next, we consider how the column density distribution changes when the same density 
field is used but all the peculiar velocities are put to zero. This is shown in Fig. ^. 

The mean transmission of the spectra computed with zero peculiar velocities differs 
from the mean transmission of the full spectra by less than a percent. It is used as the 
transmission threshold in the line identification procedure for both analyzes. The resulting 
column density distributions are very similar. 

Hence, the peculiar velocity plays a relatively minor role in determining the column 
density distribution. It changes the shapes of absorption troughs without altering the overall 
number of lines and their column densities. This serves to motivate an approximation we 
will introduce in the next section. 

A final note on velocity caustics. The reader might worry that at a velocity caustic, 
the optical depth may diverge while it is clear from equation (|T]) that for a finite number 
density of neutral hydrogen, the optical depth should always be a finite quantity. The 
resolution is that close to a velocity caustic at m = Uc, du/dx goes like \u — ud^^"^ (provided 
the second derivative of u with respect to x is nonzero, otherwise it will be (u — UcY^^ if the 
third derivative does not vanish and so on, by simple Taylor series expansion; see Shandarin 
and Zel'dovich 1989 for a similar analysis apphed to real caustics as opposed to velocity 
caustics; orbit-crossing occurs in the former but not the later). So under the integration in 
equation (^, the optical depth remains finite. We note also that because of the singular 
nature of (du/dx)'^ around u = Uc, the derivation leading to equation (|^) breaks down at a 
velocity caustic. 
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5. The Statistics of Density Peaks 



In this section we explore a simple approximation in which each density peak in x space 
is identified with an absorption line. This is motivated by the fact that peculiar velocities 
do not play a major role in determining the column density distribution and that each 
maximum in density corresponds to a minimum in transmission or maximum in absorption 
if the peculiar velocities are set to zero and if the maximum in density is separated from 
other maxima by a distance larger than that given by the thermal broadening width. 



To calculate (fNi^ya/dNm/dz, we relate dz and dx by ignoring peculiar velocities: 



dz = c ^Hdx. Hence 



Lya 



d^N, 



pk 



(24) 



dNuidz dNuidx H 

where dN^]^/ dN-m/ dx is the average comoving number density of density peaks along a 
random line of sight per unit column density, H is the Hubble constant at the redshift of 
interest . 

For each density peak, we need a simple prescription for assigning a column density. 
To that end, we perform the following expansion around each density maximum: 



ln[raHi(a;)] = ln[nHi(a;pk)] + 



1 rf^ln[nHi] 



dx"^ 



[X - Xpkj 



(25) 



It is a straightforward Taylor expansion around the position of the peak Xpk. The second 
derivative in the last term is negative. The rationale behind expanding ln[nHi] rather 
than nni itself is that nni is supposed to fall off quickly far away from the peak (until, of 
course, another peak is encountered). In other words, the above expansion implies that 
riHi has a Gaussian fall-off (instead of a power-law one if nni itself were Taylor expanded). 
In a sense, this is close in spirit to the Voigt-profile-fitting-technique. Suppose the 
broad-maximum-limit (eq. 0) holds so that the local optical depth is simply proportional 
to the number density of neutral hydrogen if one ignores peculiar velocities. Then, fitting a 
minimum in optical depth with the Voigt or thermal profile (eq. 0) is equivalent to fitting 
the corresponding neutral hydrogen density peak with a Gaussian. 



We then assign the following column density to the density peak: 



HI 



dx 

pkl + z 



nm{x) = nHi(a;pk)(l + z) ^ 



\ 



2tt 



In 



dx'^ 



["-hi] \ 



(26) 



where equation (^) has been used and where /pj^ denotes integration around the peak 
until it decays sufficiently. All x dependent terms on the right hand side are evaluated at 
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X = Xpk- The above equation is also derived by the authors in a separate paper, using the 
Stationary Phase Method (Gnedin and Hui 1996). 



Using equations ([l7|) and ([l9l), the above can be rewritten as: 



HI 



1.63 X lO^^cm-2 
/2 - 0.7(7 - 1)^ 



1.65 



-0.5 

A 



-0.7 



0.0125 



(15 



1 + z 



^\ 5 



(27) 



where A is defined by 



A= (l + 5b)'-°-'^^-') 



-t/^ In [1 + ,5b] 



(2^ 



with X being measured in Mpc. 



We will refer to our method as the Density-Peak-Ansatz. It consists of two parts: 1. 
associate each density peak in x space with an absorption line; 2. assign a column density 
to each density peak according to equation (P7|).p^ 

Making use of equations (p^) and (pT]), the column density distribution can be written 



as 



Lya 



dNuidz 



6.25 X 10"^^cm2 ' ° 



0.7 



0.0125 



1 + Z 



(29) 



2-0-7(7-1) 
1.65 



0.5 



H dAdx ' 



where x and cH~^ are in the same unit of Mpc. Most of the factors above come from the 
scaling between A and A^hi? and c/H provides the conversion from comoving coordinate x 
to redshift z. The last factor d'^Np]^/dA/dx involves is the number density of peaks in x 
space having the quantity A within the range dA. 

Let us define ^ = ln[l + 6^]- Suppose one is given P{$,,$,',$,")d^d^'d^" which is the 
probability that ^ and its first and second derivatives with respect to x fall in the specified 
ranges at a point. Then, 



dx 



/oo fU 
du dermic = o,n 
-oo J —oo 



(30) 



Strictly speaking, care should be taken not to count peaks that are separated in velocity space by distance 
much less than the thermal width as contributing to more than one absorption line. We will address this 
later in the section. 
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where is the integral of -^^j^ over all A (Bardeen et al. 1986). 



By a change of variable and a differentiation and making use of equation (l2^), one can 
obtain 

J — oo 



dAdx [2 - 0.7(7 - 1)]^ 
where ^ should be expressed in terms of C," and A using equation (^) . 

Note that the above two equations are completely general and no assumption about the 
Gaussianity of the underlying fields has been made. The hard part is of course to come up 
with the probability function P. The one point probability distribution of just ^ or density 
has been calculated for the Zel'dovich approximation ([Kofman et al. 1994|) . We find the one 
point joint-probability distribution of density and its first and second derivatives along a 
line of sight difficult to calculate analytically for the Zel'dovich approximation. A numerical 
approach is adopted in this paper and the number of peaks is counted along random hues of 
sight in actual three-dimensional realizations. In a separate paper, we discuss an analytical 
calculation based upon not the Zel'dovich approximation but the lognormal approximation, 
where ^ is assumed to be a Gaussian random field (Gnedin and Hui 1996). A comparison 
between the two will be made here. 



5.1. Testing the Density-Peak- Ansatz 

We test the Density- Peak- Ansatz in two different ways. First, we make a scatter plot of 
the column density obtained using the Threshold-Deblending- Algorithm versus the column 
density obtained by searching for the maximum density peak that contributes to each 



absorption line identified using the threshold method and then applying equation (pTj) . The 
result is shown in Fig. ^. It shows that while the agreement is far from perfect, the column 
densities assigned using the Density-Peak- Ansatz and using the Threshold-Algorithm are 
broadly consistent. 

The important question, however, is whether the Density-Peak-Ansatz, coupled with 
the Zel'dovich approximation with appropriate initial smoothing, gives the correct number 
of absorption lines as a function of column density. We compare the column density 
distribution obtained using our approximate methods against that obtained by applying 
the Voigt-profile-fitting-technique to synthetic spectra from a full hydrodynamic simulation 
(see discussion below and Zhang et al. 1996 for details). This is done in Fig. |^. Note 
that the values of 7 and Tq given in the caption of Fig. ^ is obtained directly from the 
hydrodynamic simulation. The temperature-density relation is not an exact power law but 
is well approximated by one for the relevant range of densities. We also put in the same 
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figure the column density distribution obtained using the Threshold-Algorithm, coupled 
with the Zel'dovich approximation. The predictions of the lognormal approximation are 
shown as well for comparison. 

The level of agreement between the exact hydrodynamic computation and our 
calculation based on the Density-Peak-Ansatz coupled with the Zel'dovich approximation 
is encouraging. Two sets of points are shown for our approximate calculation using the 
Density-Peak-Ansatz, one (open triangles) with exactly the same box size and grid spacing 
as the hydrodynamic simulation and the other (open squares) with larger box size and 
smaller grid spacing. They both agree very well with the exact computation. We will 
explore the effects of changing the resolution in the next sub-section. A third set of points 
(crosses) shows that the Threshold-Algorithm described in Sec. underestimates the 
number of lines at low column densities. 

The agreement between the results of the hydrodynamic simulation and our 
approximation method is perhaps telling us something interesting about the low column 
density systems. Note that for the Zel'dovich approximation computations presented 
in Fig. ^, an initial smoothing scale {ks = 2.3Mpc~^) larger than the Jeans length 
[kj = 8.4Mpc~^) was chosen to deal with the problem of orbit- crossing. The agreement can 
then be understood as a result of relatively little contribution to the low column density 
systems (10^^'^ — 10^^'^ cm~^) from structure on small length scales suppressed by our fcg ^ 
smoothing even though those scales could be larger than the Jeans length. The low column 
density systems mostly consist of relatively broad density peaks. 

Some explanation is also in order regarding the spectral analysis method used by 
Zhang et al. (1996) in obtaining the distribution shown in Fig. |^. It consists of identifying 
absorption features above a specified opacity and deblending them into individual hues 
centered at local maxima in optical depth and fitting each with a suitable Voigt profile. 
This procedure is designed to be similar to most observers' analysis techniques but is not 
exactly the same. In particular most observers' method consists of seeking a "minimal" set 
of Voigt profiles, the superposition of which reproduces the input absorption spectrum to 
within some specified error consistent with noise. It should be noted, however, that the 
definition of "minimal" tends to vary from one observer to another. The superposition of 
Voigt profiles identified using the method of Zhang et al. seems to reproduce the input 
spectrum quite well (see Fig. 2 of Zhang et al. 1996), although the agreement has not been 
carefully quantified. One should keep in mind possible differences in the column density 
distributions obtained using different methods (see Dave et al. 1997 for related discussions). 

We also show in Fig. ^ two sets of curves based on the lognormal approximation but 
using the same Density-Peak-Ansatz (see Gnedin and Hui 1996). One of them has the 
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same amount of initial smoothing as that of the truncated Zel'dovich approximation and 
the other has less smoothing so as to match the final (not linear) rms density fluctuation of 
the Zel'dovich computation. In both cases, the lognormal approximation tends to predict 
too much flattening of the column density distribution at low column densities. In general 
the lognormal approximation tends to give column density distributions that deviate quite 
significantly from power-law unless a very small smoothing length k^^ is chosen. (From 
the figure, it might appear that the lognormal approximation gives more lines than the 
Zel'dovich approximation at the very low column densities but it is really a resolution effect: 
see the next sub-section.) 

The reader might have noticed that we have included in Fig. ^ a wider range of column 
densities than is warranted by the nature of our approximations. For instance, objects with 
column densities higher than 10^^ cm~"^ are almost certainly highly nonlinear and we do 
not expect the truncated Zel'dovich approximation to work well in this regime. For the 
low column densities, the finite resolution should cause us to underestimate the number 
of absorption lines. In the next sub-section, we give quantitative estimates of the range 
of column densities within which the Density-Peak-Ansatz, used in conjunction with the 
truncated Zel'dovich approximation, can be counted upon to give reliable column density 
distributions. 

We will also discuss two different ways of defining a density peak in the next sub-section: 
local maxima over three cells or local maxima over five cells with the slope on either side 
of the maxima not changing signs. The three-cell criterion is used in Fig. ^. One expects 
however that some of the three-cell peaks are not real but merely artifacts of numerical 
noise, especially those with low column densities. The five-cell criterion, on the other hand, 
probably fails to include some narrow peaks which are real. We will see that the two 
different criteria give almost identical results above a certain column density. 

One aspect of the Density-Peak-Ansatz we have glossed over is that two density peaks 
that are separated by a distance in velocity space much less than the thermal width should 
be counted as contributing to not two but one absorption line. A more sophisticated 
approach would be to group together such density peaks and use the sum of their column 
densities as the column density of one single absorption line. We find that for the the range 
of validity discussed in the following sub-section, it makes little difference. It is conceivable, 
however, that this effect cannot be ignored for simulations with higher resolution than we 
have, or at higher redshifts, where line blending is more important. 
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5.2. The Range of Validity 

For the computation presented in Fig. the column density (given by the Density- 
Peak- Ansatz) above which the mean exceeds 5 is about 10^^'^ cm~^. For the parameters 
listed in the caption of Fig. |, A^hi = 3.6 x 10^^ A cm~^ (eq. |^). We therefore take 
A = 3.5 as an upper limit beyond which we cannot expect our approximations to be reliable. 
Shock-heating should be relatively unimportant for (5b less than about 5 (see Hui & Gnedin 
1996). 

Note that according to Fig. ^, comparing with the hydrodynamic simulation data, 
the Density-Peak-Ansatz coupled with the truncated Zel'dovich approximation, seems to 
give a reliable number density of absorption lines even for column densities higher than 
lO^^'^cm"^. The level of agreement at such high column densities (and by extension, such 
high 5b) is surprising. We will adopt the conservative upper limit of A = 3.5. 

To determine the column density below which finite resolution results in an 
underestimate of the number of absorption lines, we perform a simulation using the 
truncated Zel'dovich approximation with the same parameters as the open squares in Fig. ^ 
but higher resolution: same box size of 12.8 Mpc but smaller grid spacing of 0.0284 Mpc. 
A comparison of the resulting column density distributions is made in Fig. ^. 

Note that we have included two definitions of density peaks (three-cell and five-cell). 
For each simulation, the true column density distribution is probably somewhere between 
the two in the places they differ. 

We take the low column density cut-off to be 10^^'^ cm^^ for the lower resolution 
simulation (box size of 12.8 Mpc , with grid spacing of 0.05 Mpc) using the three-cell 
definition of peaks. It can be seen that the higher resolution simulation differs from the 
lower one only at column densities less than roughly this cut-off value. Moreover, above 
this column density, the three-cell and five-cell criteria give almost identical results. 

The parameters in the simulations in Fig. |^ are such that A^hi = 3.6 x lO^^Acm"^ (eq. 
||27||). Hence the above column density cut-off implies a lower limit of 0.18 for A. From now 
on, we will use the three-cell definition of density peaks. 

For readers interested in applying our formalism at different redshifts, we recommend 
that they choose the appropriate range of A using the same methods above: upper limit 
set by non-linearity and lower limit set by running simulations of varying resolutions. The 
range is not expected to change significantly with redshift. 

In the following section, we systematically investigate how the column density 
distribution depends on the cosmological parameters and properties of the intergalactic 
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medium. All the simulations presented in the next two sections have the same resolution 
and box size, 256^ grid points with grid spacing of 0.05 Mpc. For each of them, we 
will only plot the part of the column density distribution that falls within the limits of 
0.18 < A < 3.5. The column densities these limits correspond to depend on the properties 
of the intergalactic medium and the redshift (eq. Note that our conservative limits for 

A greatly reduce the range of column densities we can examine, but within these limits we 
can be reasonably confident that the truncated Zel'dovich approximation together with the 
Density-Peak-Ansatz should yield accurate predictions for the column density distribution. 



6. The Column Density Distribution: Dependence on the Ionization Flux, 
Temperature, Equation of State and the Mean Baryon Density 

We use the Cold Dark Matter (CDM) model presented in Fig. ^ to study systematically 
how the column density distribution depends on properties of the intergalactic medium: 
the level of radiation background, the mean baryon density and its thermal properties. The 
tools we use to calculate the column density distributions are the truncated Zel'dovich 
approximation and the Density-Peak-Ansatz. 

The transfer functions for this model and all the other models considered in this paper 
are taken from Ma (1996). A summary of parameters of all models can be found in Tables |I| 
and H 



6.1. Dependence on Overall Temperature, Ionization Flux and Baryon Density 



Let us first consider the case in which the equation of state is fixed at T oc (1 + (5b)°'^. 
As is shown in equation (pT]), the column density of a density peak with a given 5b 
(overdensity) is proportional to the following combination of parameters: 

F = ( ] I 1 (^)' ■ (32) 



0.0125 



Hence, by equation (|30|) , if F is rescaled by a certain factor (by changing Tq, f2b or Jhi or 
their combinations), the number of absorption lines is also changed by the same factor at 
an appropriately rescaled column density. 



Suppose F is rescaled to F' such that F' 

N' 



dN^jidz 



rF, then 

1 C^^iVLya 



r dNmdz 



(33) 
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It implies that if the column density distribution is a pure power law, then in a log-log 
plot of the number of absorption lines per unit column density per unit redshift versus 
column density, the straight line would simply be shifted to the right or left (or up/down) 
by rescaling. In reality, the column density distribution only approximately obeys a power 
law and so there should be a slight change of slope at any given column density as a result 
of rescaling. 

The effects of rescaling can be seen clearly in Fig. |^, where F is allowed to take the 
values 0.25, 1 and 5. Keeping VL^h'^ = 0.0125 and Jhi = 0.5, it corresponds to changing Tq 
from about 72000 K to 1000 K. Alternatively, keeping Tq and Qbh"^ fixed at their canonical 
values (as shown in eq. |^), it corresponds to allowing Jhi to vary between 2 and 0.1. See 
Hui and Gnedin (1996) for a discussion of the dependence of Tq on reionization history. Tq 
is expected to fall within the range quoted above. 

The conventional value of Qbh"^ = 0.0125 has been challenged by recent measurements 
of light element abundance in high redshift absorption systems. Tytler and Buries (1996) 
obtain a value of 0.024, which for Tq = 10^ K and Jhi = 0.5, implies F = 3.686, well within 
the range of F plotted in Fig. ^ The analysis of Rugers and Hogan (1996), on the other 
hand, favors the value 0.006, which means F = 0.23 for the same values of T and Jhi. The 
lowest set of points in Fig. ^ has to be lowered further to accommodate this value of the 
baryon density. 

The observational data are taken from Hu et al. (1995), measured at about redshift of 3 
and corrected for incompleteness. The incompleteness correction was obtained by applying 
the same spectral analysis technique for the observed data to simulations of randomly 
distributed Voigt-profiles with a known column density distribution and measuring how 
much the analysis method underestimates the number of low column density lines. The 
amount ranges from about no correction necessary at A^hi > 2 x 10^^ cm~^ to a factor of 5 
increase in the number of lines at N^j ~ 3 x 10^^ cm^^. 

We note in passing that strictly speaking, altering Qh, in addition to rescaling the 
number of absorption lines as discussed above, also changes the transfer function in 
a non-trivial way. But the effect is very small for models in which the dark matter 
(non-baryons) dominate. In fact, using the BBKS ( Bardeen et al. 1986 ) transfer function. 



which does not take into account the effect of baryons at all, instead of the Ma (1996) 
transfer function, which does take it into account for Q^h"^ = 0.0125 with h = 0.5, gives 
essentially the same column density distribution for the range of column densities discussed 
here. For models where the baryon content is proportionally higher. Low-density Cold Dark 
Matter Models for instance, changing Qbh"^ has a more pronounced effect on the transfer 
function. 
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6.2. Dependence on the Equation of State or the Temperature-density 

Relation 

Let us hold fixed Tq, flhh^ and Jhi at their canonical values as shown in equation 
(^) but allow the equation of state to change, for the same CDM model as above. As 
is pointed out in Sec. |3.2| , the temperature-density relation for low enough overdensity is 
well-approximated by a power law where the power index is around 0.5, but can change 
slightly depending on the reionization history. We plot in Fig. ^ the column density 
distributions for 7 = 1.2, 1.5, 1.7 where 7 is defined by T oc (1 + 6^)^'^. It should adequately 
cover the possible range of 7 ( [Hui fc Gnedin 19961 ). 

The first thing to notice is that the column density distribution remains almost the 
same for the three different values of 7. This is because 7 affects column density through 
the power index of (1 + 5b), which is 2 — 0.7(7 — 1) (^q- EH)- The index does not change 
significantly for the range of 7 considered. A larger index (smaller 7) means for a density 
peak with a given 1 + 5b (and its second derivative), the column density is larger or smaller 
depending on whether 1 + (5b is bigger or smaller than one. The net effect is to decrease 
the slope of the column density distribution. The effect, though very small for the values 
of 7 plotted, can still be seen in Fig. |. We also show the approximate slopes given by an 
analytical formula (eq. ^) which will be discussed later. Note how the column density 
distribution does not exactly follow a power law but can be approximated by one. 

Hence as a crude approximation, we conclude that the mean temperature, radiation 
intensity and baryon density mainly determine the overall normalization of the column 
density distribution. The equation of state, on the other hand, mostly affects the slope of 
the column density distribution, but its effect is small for a reasonable range of 7. 

7. The Slope of the Column Density Distribution 

It has been shown that the normalization of the column density distribution is 
influenced by the thermal and ionization states of the intergalactic medium, which are not 
well-constrained observationally. The slope of the distribution, on the other hand, is only 
weakly affected by the equation of state or temperature-density relation. 

The slope of the column density distribution is therefore relatively free of uncertainties 
due to our ignorance of the thermal and ionization properties of the intergalactic medium. 
We now turn our attention to the effect of the power spectrum on the slope of the 
distribution. 
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From equations ( ]27|) and (pS]), it can be seen that the column density A'hi is proportional 
to (1 + Sby~'^''^^"'~^^ times which basically defines a length scale. Taking into account 

the correlation between this length scale and the overdensity, we find from our simulations 
(which use the Zel'dovich approximation) a useful approximate relation for column densities 
between about 10^^'^ and 10^^'^ cm~^ (or more accurately, for the range of validity discussed 
in Sec. W^ '- 

iVHioc(l + 5,)i-^«-°-^(^-i), (34) 
which roughly means the length scale 1/a/^ is approximately proportional Q to {l + 6b)^^'^'^ ■ 

Since we are interested in the slope of the column density distribution, the relevant 
quantity to consider is: 

Ti ^_ d\nf_^dew'm,c = o,n .... 

dlnNni ~ ^1.68-0.7(7-1)' d^ ' ^^^^ 



The equality follows from equations (|3^), (p^), (|30D and (|3lD and noting that ^ = 1 + 5b. 
The column density distribution can be approximated by the simple power law N^f if (3 
defined above is only weakly dependent on C, or N-m. 

Lacking an analytical expression for P under the Zel'dovich approximation, we can 
nonetheless guess what the general properties of the quantity m are. First of all, m depends 
on ^ in general because the integral d$,"\^"\P{^, ^' = 0, ^") cannot be a simple power law 
in ^. This is because we expect the integral to vanish for very large and very small ^ and 
peak at some intermediate ^. This implies one should not expect an exact power-law for 
the column density distribution, although pieces of it might be approximated by power-law. 
Suppose ,^pk is the value of ^ where the integral jloc,d^"\^"\-Pi^^^' — O^'C") reaches its 
maximum value, then the column densities we are interested in must correspond to ^ > 
where ciln jloo'^^"\^"\-P{^^^' — 0;^")/^^ is negative. This is based on the knowledge that 
the computed (as well as observed) slope in equation (^) is less than —1 (the factor 
1.68 — 0.7[7 — 1] is always positive for reasonable values of 7). 

Under the Zel'dovich approximation, the quantities ^ and its derivatives can be related 
to the displacement potential ipi^q) using equations (|TT|) and ( [T^ ) and so one can express 
P(^, C,'.^") in terms of the probability density for derivatives of ip: 

f de\i"\P{^, = 0, = C <#.,#i,fc#.,fcde"|e-^P(^,„ ij^.ui) (36) 
J —00 J 00 



^^A log-log plot of 1/%/^ versus (1 + Sb) actually shows a lot of scatter but equation ( |3^ ) appears to 
capture the overall dependence of A^hi on 1 + St. 
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where Sd denotes the Dirac delta function, ipij is the derivative of the displacement potential 
ip with respect to and qj (eq. |]1T|) and similarly for ipij^ and ipijki and e* is the unit vector 



pointing along the line sight. The quantity P{ipij, 'ipijk, i^ijki) is the probability density of ipij 
and its derivatives in Lagrangian space and the factor converts it into its counterpart 
in Eulerian space (see Kofman et al. 1994). The quantities ^ and its derivatives along the 
line of sight can be expressed as functions of ipij, ipijk, ipijki and e\ Isotropy of the universe 
implies that one can average the above expression over all possible orientations of e*. 

For the cosmological models we consider in this paper, the probability density 
P{4'ijy4'ijk,i'ijki) is a multivariate Gaussian function that depends on 3 parameters: ctq, <Ji 
and (J2 defined as follows (see Bardeen et al. 1986): 



D+{t)J A7rk^+^^P{k)e-W''srdk, j= 0,1,2 (37) 



where is the linear growth factor which is equal to {1 + z)~^ for a universe at critical 
matter density and ks is an appropriate smoothing scale. This follows from the structure of 
the various expectation values: < ipijipmn >c<: Cq, < ipijipmnip >(x erf, < ipijki'ipmnop >oc cr| 
and < tpiju-ipmni >oc af. The fact that < ipiji^mni > and < ipijki-ipmno > vanish by 
isotropy means that the probability density factors into two separate multivariate Gaussian 
functions: P{^pij,iJijk,ipijki) = P{ipij,i^ijki)P{'iptjk), with the first factor depending on ao, ai 
and (72 and the second factor depending on (Xi. 

One can replace ai and a2 by combinations of ctq and the following two new parameters, 
which Bardeen et al. (1986) defined: 

2 

= 73— , = ^ , (38) 
CT2 a2ao 

where we have renamed 7^ to distinguish it from 7 we use in this paper. 

While R^, defines a length scale, the quantity 7^ is a measure of the slope of the power 
spectrum. We find it convenient to use the following quantity rzefj in place of ■jb (Bardeen 
et al. 1986): 

57I-3 

It is easy to show that the above quantity coincides exactly with the slope of the power 
spectrum if it obeys a pure power-law. For 7^ = 0.5, rics = —2.33. The riefr defined here 
should be distinguished from n in Tables |I] and H: n is the slope of the power spectrum at 
large scales whereas Ues is the slope at the (small) smoothing scale. 

The advantage of the new notation is that it is possible to show, by changing variables 
ipijk = ipijkR*, ipijki = i^ijkiRl and ^" = C,"Rl, the integrals in equation (|37|) are independent 



nefT = . (39) 
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of except for a normalization factor. This implies, by equations (^) and (|37|), that 
the slope {3 depends on the equation of state through 7 and on the power spectrum only 
through (Jo and ti^q. 

It is hard to make analytic progress from this point on, because ^ and its derivatives 
are complicated functions of the derivatives of if). We resort to our numerical simulations to 
extract the parameter dependence of (3. Using the above arguments and observing that for 
all models we consider in this paper, ctq ~ 1 and n^Q ~ —2.33 (7^ ~ 0.5), we assume the 
following form for (3: 

gi + i?2(ao-l) + i?3Kff + 2.33) 
^ 1.68-0.7(7-1) ' ^ ' 

where B2 and are constant coefficients of a Taylor series expansion of m (eq. 

We determine these constants by computing j3 for a series of CDM models of varying 
(To and 7b- The column density distributions are obtained by using the Zel'dovich 
approximation together with the Density-Peak- Ansatz, as described in previous sections. 
We pick out particularly models with similar rzefr's but very different a^s and vice versa. 
It is found that the following expression for (3 fits reasonably well the slope of the column 
density distributions for the CDM models, as well as the Cold + Hot Dark Matter (CHDM) 
models that we will also study later: 

0.96 - 0.8(^0 -l)-0.4Kff + 2.33) 
^ ^ 1.68-0.7(7-1) ' ^ ' 

where /3, 7, and n^a are defined in equations (|35|) , ([T9|) , (|37|) and (^) respectively. 



We show two sets of examples in Fig. ^ and Fig. The former has two CDM models 
with the same ries but different (Tq's while the opposite is true for the latter. To find models 
with such properties, we find it necessary sometimes to smooth on scales larger than the 
orbit-crossing scale k'^^ defined by equation ( |15|) (examples are models CDM3 and CDM4 
shown in the Figures; see also Table p. For this reason, the column density distributions 
computed in this section should be viewed as predictions of the corresponding CDM models 
at the specified smoothing scales only: they do not necessarily coincide with the predictions 
of these models if more appropriate smoothing scales are chosen. Nor should the CDM 
models shown in the two figures be considered realistic models of the universe. 

Examples of how well equation ( ^T] ) describes the variation with the equation of state 
for a CDM model can be found in Fig. ^. More examples for CHDM models are also shown 
in Figures |15|, [16], and |18[ 



Before we go on to discuss physical interpretations of the above expression for /3, let 
us make some general remarks and state a few caveats. First, as we discussed before, j3 
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is expected to vary with column density but can be approximated as a constant over a 
restricted range. That a power-law column density distribution is only an approximation 
is clear even in some of the figures mentioned above. We find it appropriate to allow for a 
maximum error-bar of ±0.1 for 13 in equation PTI), given 7, o"o and 7^. The coefficients -Bi, 



B2 and i?3 (eq. ^y]) should be viewed as correspondingly uncertain. Furthermore it should 
be kept in mind that equation ( PD holds only for (Tq ~ 1 and rieff ~ —2.33. as the constants 
B2 and i?3 are only meant to be coefficients of a Taylor series expansion. 

Finally, one obvious problem of using equation (^) to make prediction for a given 
cosmological model is that both cxo and neg varies with the smoothing scale. An example 
of how the choice of smoothing scale can affect the column density distribution is shown in 



Fig. O. The comparison between results of a hydrodynamic simulation and our Zel'dovich 
computation in Fig. ^ seems to support the choice of smoothing in equation (|T3p . However, 
one could reasonably argue that the Jeans scale (eq. |T^) is a more physically motivated 
smoothing scale while the choice in equation ([T5| ) was made only as a device to counter the 
effect of orbit-crossing and hence one cannot be sure that one is not thereby erasing real 
structures on scales smaller than the orbit-crossing scale but larger than the Jeans length, 
which can contribute significantly to the low column density Lyman- alpha forest. One 
way to settle this question is to make more detailed comparisons with more hydrodynamic 
simulations, which is outside the scope of the present work and we will leave for a future 
paper. It suffices to say that, for models with sufficiently little small scale power such as 
the CHDM models considered in Sec. |^, which have orbit-crossing scales close to or smaller 
than the Jeans length, it is probably safe to smooth according to the prescription laid down 



in Sec. |3.1| and that the precise smoothing scale might not matter very much (see Fig. 0). 
For more nonlinear models, we can only offer the comparison in Fig. |^ as one piece of 
evidence that smoothing on the orbit-crossing scale works reasonably well. However, for 
whatever smoothing scale one decides upon, we expect equation (|¥lD to hold approximately, 
provided that do ~ 1 and ~ —2.33. 

With these being said, let us try to understand qualitatively the parameter dependence 
of (3 as expressed in equation (|4l|) . 



First, the dependence on 7. Recall that A^m oc (1 + ^j^) 1-68-0.7(7-1) (taking into account 
the correlation between 1 + 5b and its second derivative). This means that higher density 
peaks translate into higher Ahi and vice versa. Lowering 7 increases the effect of this 
translation i.e. a given density peak with 1 + 5b > 1 is associated with a higher column 
density if 7 is reduced (the opposite is true for 1 + 5\, < 1). The net effect is to stretch a 
power-law-like column density distribution and make it flatter. The effect is small for a 
reasonable range of 7. It should be between 1.3 and 1.62 at z = 3 if the universe reionizes 
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before 2; ~ 5 (see Hui and Gnedin 1996). A somewhat larger range is shown in Fig. ^. 

Second, the dependence on cxo: a cosmological model with lower ctq is in a more "linear" 
state of evolution compared to higher ctq models. In other words, a lower (Tq model has 
proportionally more intermediate density peaks compared to high density ones, hence the 
steeper column density distribution (associating once again high density peaks with high 
column densities using the A^hi oc (1 + 5^)i-68-o-7(7-i) scaling). For sufficiently low column 
densities, however, the absorption lines arise from very underdense regions which should be 
more common in higher ctq models. Hence at very low column densities, the high uo model 
should win: it has more very low density peaks. Where this might occur we cannot tell 
from our simulations because of the limited resolution. For the range of column densities 
we can measure reliably, the slope of the column density distribution simply steepens as ao 
is lowered. 

The dependence of (3 on nefr is more subtle. To understand it, we resort to a 
Press-Schechter-type argument (Press & Schechter 1974). According to the Press-Schechter 
theory, the low mass slope of nc{M) (the number density of clumps of mass M per unit 
range dM) is — 1.5 + neff/6 which means that smaller n^^ {ries < — l)implies a steeper 
slope of the clump number density distribution, exhibiting the same trend as in the case of 
equation (pj). The n^^ dependence arises from the fact that (jQ{k) oc A;("<:ft+3)/2 (^^^ with 
smoothing scale k) and assuming M oc , together with Mnc{M) oc k^/ao{k) (for high 
k). The factor of k^ appears because nc{M) is a number density. The dependence of nc{M) 
on ao{k) follows from a simple conjecture: the probability that a given point belongs to a 
clump of mass M or higher is equal to 2 P{6)d6 where 6th is a fixed threshold and P{S) 
is the Gaussian probability distribution for overdensity 6 with dispersion o"o(fc). 

One can carry the above reasoning over to absoprtion lines, keeping in mind something 
like the Threshold Algorithm for identifying lines as the analogue of the threshold-criterion 
of the Press-Schechter theory for a collapsed clump. Assume A^hi oc k~^ (because in the 
linear regime, it is the peak width, not the peak height that provides the lowest order 
contribution to column density), it "follows" that Nmi^d"^ N^ya/ dN^dz) oc k/ao{k) , again 
in the high k limit. The factor of k is to account for the quantity of interest being a 
one-dimensional number density. It is easy to see this implies /3 = (1 — nefr)/2 i.e. smaller 
rieff (larger |nefr|) implies steeper column density distribution. This admittedly crude 
argument actually gives a value for (3 that agrees surprisingly well with what we have found 
using the Zel'dovich approximation and Density-Peak-Ansatz. For example, for the CDM 
model in Fig. ^ rics = —2 (see Table |ip, P = 1.5 according to the adapted-Press-Schechter 
argument. This of course brings up the question why there does not seem to be any 
(To dependence of /? in the above argument. It is in fact possible to develop this line of 
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reasoning further to include the cxo dependence, but we will leave this to a later paper. 

Finally, it is interesting to point out that because most realistic cosmological models 
have similar power on large scales (COBE or cluster scale), small (Jq on small scales almost 
necessarily means a steeper fall-off of the power spectrum, namely smaller Ues- So, in 
practice, both the amplitude cxo and the slope Ues of the power spectrum work in the same 
direction in decreasing/increasing the slope of the column density distribution. Hence, as a 
rule of thumb, less small scale power implies a steeper column density distribution, for the 
column densities of interest in this paper. 



8. The Column Density Distribution for CHDM models 

In this section, we turn our attention to CHDM models, which have relatively little 
small scale power and for which the method of truncated Zel'dovich approximation is 
ideally suited. The orbit-crossing scales for them are either smaller than or only slightly 
larger than the Jeans length and so one has reasons to believe any structure erased by our 
smoothing procedure does not contribute significantly to the number of absorption lines in 
the range of column densities we are interested in. Predictions for the CHDM models are 
also of current interest because no hydrodynamic simulations of the Lyman-alpha Forest 
have been performed for these models. 

We list in Table ^ the CHDM models considered in this section. They are all f2o = 1 
models with = 0.05. Both the = 0.2 and = 0.1 versions are considered. They 
have been shown to give good agreement with the observational data on large scales {k 
around 0.02 — 0.4 Mpc~^) (see Fig. 6 and 7 of Ma 1996). The Qi, = 0.3 models seem to 
conflict with the observed abundance of damped Lya systems, which correspond to roughly 
k around 0.1 — l.OMpc"^ comoving in the linear power spectrum (|Mo fc Miralda-Escudej 
19941 ; IKauttmann fc Chariot l994| ; |Ma fc Bertschinger 199^ ). We include one fi^ = 0.3 



CHDM model in Table ^ for the sake of comparison. As is shown convincingly by Ma 
(1996), all models need some amount of tilt to match observations. 

We compute as before the column density distribution for each model using the 
Density-Peak-Ansatz and the Zel'dovich approximation with appropriate smoothing. The 
(density weighted) power spectrum for each CHDM model is taken from Ma (1996). 

As we have discussed in the last section, a plot of rms smoothed linear density 
fluctuation do versus smoothing scale ks (eq. [^) is a very good indicator of what column 
density distribution to expect. This is done in Fig. |T^ for the CHDM models tabulated in 
Table 0. The no-tilt as = 0.7 CDM model is also plotted for comparison. 
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Because of neutrino free streaming, all CHDM models have less power (smaller (Tq) and 
steeper spectral slope (smaller rics) than the CDM model on small scales. Those with more 
neutrino content {Qi, = 0.2) have even less power than the others. In fact, the = 0.2 
models have cxo < 1 on all scales larger than the Jeans scale (fcs < fcj ~ lOMpc"^). One 
expects the Zel'dovich approximation to work particularly well for these models because 
the amount of orbit-crossing will not be significant, even without initial truncation. 

This is borne out by the next test: we compute the column density distribution for 
one CHDM model B2 and examine the effect of choosing different smoothing scales. The 



result is plotted in Fig. |T3|. The column density distribution in the range plotted does not 



change much at all for the three different smoothing scales plotted. Contrast this with the 



case of (Js = 0.7 CDM (Fig. |ri|) where the column density distribution is more sensitive to 
changes in the smoothing scale. That's why the truncation scale has to be chosen with some 
care for the CDM model: not too small [k^ too big) so that too much orbit crossing has 
occurred and not too large so that too much small scale structure is erased. The standard 
prescription (Sec. jSTTI) seems to work well according to Fig. |^. 

For the CHDM model considered (in fact, it holds true for all other Q,y = 0.2 models 
here), the amount of small scale power is so insignificant that excluding them by smoothing 
does not affect the overall column density distribution at all (except possibly that one loses 
the small scale fluctuations that can give rise to very low column density absorption i.e. 
lower than our resolution limit). We have also done similar tests for the fij, = 0.1 models, 
their response to changes in the truncation scale is somewhere between the as = 0.7 CDM 
model and the fli^ = 0.2 CHDM models, as can be expected based on their difference in 
Fig.O. 



We adopt the following truncation scales for the CHDM models. For the Q^, = 0.1 
models, the standard prescription described in Sec. ^]T| is used (i.e. ks = 1.5 /cnl)- The 
flu = 0.2 models, according to the above prescription, would have truncation scales less 
than the Jeans length {ks > kj) and so by the arguments presented in Sec. [Ol , ks = kj is 
adopted. Again, we emphasize that for this class of models that have relatively little power 
on small scales, the precise truncation scale is not important. A summary of the truncation 
scales for all models can be found in Table 0. 



The CHDM models with fij, = 0.1 are plotted in Fig. IT3|. Values of F that give 



reasonable match to the observational data are chosen for each model. Note how the 



low-Hubble-constant-models {h = 0.5) require a slightly lower F (eq. ||32|) than the 
higher-Hubble-constant-models. The equation of state is chosen to be the same for all 
models (7 = 1.5, see equation [|1^). The level of agreement with the observational data. 



for the given choices of parameters, is satisfactory. Notice how the low Hubble-constant 
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{h = 0.5) models tend to have steeper column density distributions, because they have 
less power on the relevant scales (see Fig. |T2|). Their slopes can be brought into better 
agreement with that of the observational data if a smaller 7 is used. 

For the Qi, = 0.2 models, we cannot find values of F that give the same level of 
agreement with observations for 7 = 1.5. Two examples are shown in Fig. |16| and Fig. 
Both have h = 0.5 and small amounts of tilt. For each, three sets of theoretical predictions 
are plotted, one for each value of F: 1, 2.5 or 5. For Q^h"^ = 0.0125, the conventional 
Big-bang nucleosynthesis value, and To = lO'^K, they correspond to radiation intensity Jm 
of 0.5, 0.2 and 0.1 (eq. [^). As we have shown before, changing F mainly shifts the sets of 
points without altering the slope significantly. For the column density between about 10^^'^ 
and lO^^'^cm"^, the slope of the predicted distribution seems to be too steep compared to 
the observational data. (3 is about 1.86, with some flattening at the lower column densities, 
compared with the observed value of about 1.5. 

Another = 0.2 CHDM model {D2 in Table 0), which has a higher Hubble constant 
{h = 0.65), is shown in Fig. [1^. (The C2 CHDM model which also has h = 0.65, gives 
very similar column density distribution.) The slope of its column density distribution is 
not as steep as the previous ones. This is expected because the higher Hubble constant 



models have slightly more power on relevant scales, as is evident in Fig. |I2|. In fact, one 



might argue that the middle set of points in Fig. |1^, the one having F = 3.57, matches 
the observational data reasonably well if both observational and theoretically errors are 
carefully taken into account. However, it is still true these two models predict a steeper 
column density distribution for A^hi between about 10^^'^ and lO^^'^cm"^, compared to the 



n^ = 0.1 CHDM models (Fig. |T4D 



It is not hard to understand the column density distributions of the CHDM models 
presented if one refers back to Fig. or Table 0. The fli^ = 0.2 models have less power 
(smaller ctq) and steeper spectral slope (larger Inefrl) than those with fl^ = 0.1 on scales 
lMpc-^<fcs<10Mpc-\ which are relevant for the range of column densities we are interested 
in. As we have explained before, the column density distributions are therefore steeper for 
the Qi, = 0.2 models at this range of column densities. Among the = 0.2 models, those 
with a lower Hubble constant produce comparatively steeper column density distributions 
because they have even less small scale power and steeper spectral slope than the ones with 
a higher Hubble constant. 

In discussing the predictions of the slope of the column density distribution for these 
models, one should bear in mind uncertainties due to the choice of smoothing scale, as 
we have presented in the last section. However, the fact that these CHDM models have 
relatively little power on small scales works in our favor: the exact choice of the smoothing 
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scale does not affect the slope of the distribution very much (see Fig. ^). This is especially 
true for the Q^, = 0.2 models which have orbit- crossing scales (eq. |T5|) smaller than the 
Jeans length. So, the general conclusion of steeper column density distribution for the 
tilted CHDM models considered here compared to, say, the CDMla model (see Table ^ and 
Fig. 1^; the results of which compare favorably with a hydrodynamic simulation) should be 
robust. 

In Sec. we have discussed how the equation of state or temperature-density relation 
can also change the slope of the column density distribution, although the effect is small 
for a reasonable range of 7. We show in Fig. |1^ the effects of altering the equation of state 
on the column density distribution for one particular CHDM model {A2). F is fixed at 
2.5, the value that seems to give a column density distribution closest to the observational 
data. Smaller 7, as we have noted before, helps flatten the column density distribution but 
the flattening seems to be not quite enough even for 7 = 1.2. We show in the same figure 
a dashed line with a slope of —1.75 (which follows from equation (|4T|) by putting 7 = 1.2 
and substituting the appropriate values of cxo and ries as given in Table H). It should be 
kept in mind that we can always shift the column density distribution up and down by 
rescaling F (Sec. P), so the normalization is not important. It seems 7 < 1.2 is needed for 
this model to give the right slope of the distribution, at least the right slope to within the 
95% confidence limits of the observed (3 (1.37,1.51). The same conclusion holds for the 
other low Hubble constant fij, = 0.2 model {B2). We should emphasize, however, that a 
more detailed comparison between the predictions of the models and observations, taking 
into account noise and biases of the line identification techniques, is necessary before any 
model can be considered ruled out. 

The high Hubble constant {h = 0.65) Qi, = 0.2 models C2 and D2, on the other 
hand, have intrinsically flatter distributions and a reasonable match between theory and 
observations can be made by choosing 7 in the range 1.2 — 1.7. 

9. Conclusion 

We have systematically developed a set of tools to compute in an efficient manner the 
column density distribution given a cosmological model. One fundamental assumption of 
the approximations involved is that most of the Lja forest with column densities in the 
range 10^^'^ — lO^'''^ cm~^ originates from regions of low overdensities or even underdensities 
which have not undergone orbit-crossing. The result of a comparison with a hydrodynamic 
simulation lends support to it. 
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One major conclusion we reach, in the process of developing the tools, is that the 
peculiar velocities play a minor role in determining the column density distribution at 
our column densities of interest, even though they are very important in determining the 
shapes of individual absorption line profiles. We take advantage of this fact and develop 
a method we call the Density-Peak- Ansatz in which each density peak is identified as an 
absorption line and assigned a column density based on its local properties. The column 
density distribution then becomes a statistic of density peaks. 

In Sec. P and |^, we investigate the factors controlling the column density distribution, 
which can be divided into two categories. One mostly affects the normalization, while 
the other mostly influences the slope. Those that fall into the former category include 
the ionizing radiation intensity, the mean temperature of the intergalactic medium and 
the mean baryon density. Uncertainties in their values are such that almost any viable 
cosmological model which has the correct slope of the column density distribution can be 
made to match observations by a judicious choice of parameters. 

The factors that mostly affect the slope of the distribution include the equation of 
state and more strongly so, the amplitude and slope of the (linear) power spectrum on 
scales 1 Mpc~^ Hk^ lOMpc^^. Models that have less power on these scales tend to have 
comparatively more intermediate density peaks than high density ones and hence have 
relatively steeper column density distributions. Models with a steeper power spectrum on 
these scales have fc^/(To(fc) that falls off more quickly with increasing (o"o(fc) oc k^"'''fi+3)/'i^ 
and therefore a steeper column density distribution, assuming high column density 
corresponds to high (Press-Schechter argument; see Sec. Equations of state which 
are closer to isothermal (smaller 7 where 7 satisfies T oc (1 + 6\y)'^~^) tend to produce 
flatter column density distributions. However, within the reasonable range of 7 (see Hui 
and Gnedin 1996), its precise value depending upon the reionization history, the effect of 
changing the equation of state is small. We put forward an approximate expression relating 
the slope of the column density distribution to 7 and the amplitude (do) and slope (rics) 
of the power spectrum on small scales (eq. ^l|), which describes reasonably well all the 
models we study in this paper. 

Hence, the slope of the column density distribution provides a measure of the amplitude 
and slope of the power spectrum on small scales for a given cosmological model and given 
temperature-density relation. We apply our techniques to study a class of CHDM models 
which are known to have less power on small scales compared to other popular CDM models. 
We conclude that the CHDM models indeed produce steeper column density distributions 
compared to the CDM models. In particular, the low Hubble constant {h = 0.5) = 0.2 
CHDM models, which have the least amount of power on small scales among the models we 
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study, have column density distributions which can be made consistent with observations 
only for 7 less than the values we consider reasonable. We emphasize however that only 
after a more detailed comparison between theories and observations, including all the effects 
of noise and biases of the line-identification methods, can any model be considered ruled 
out by the observed column density distribution. 

We therefore conclude that a lot of work still needs to be done both on the 
observational and theoretical fronts. The biases of the line-identification techniques used 
for data reduction deserve close study so that the error bars in the observed column density 
distributions can be better understood and perhaps reduced. Numerical simulations on the 
CHDM models should be carried out to further test the accuracy of the approximations 
made in the present work. Detailed comparisons with hydrodynamic simulations will shed 
light on the appropriate choice of smoothing scales for approximate methods such as the 
one we present here. The effect of a fluctuating radiation field, instead of a uniform one as 
is assumed here, has to be investigated. Moreover, in terms of constraining models, it is 
also important to examine other possible statistics. We have shown, for instance, that the 
column density distribution is relatively independent of peculiar velocities. Are there other 
statistics that can take advantage of the different peculiar velocity structures predicted by 
different cosmological models? 

In short, the study of the Lya forest has entered an exciting stage. There is a gold 
mine of information contained in the quasar absorption spectra waiting to be discovered. 
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A. Smoothing at Jeans Scale 

The effect of gas pressure is to smooth the baryon density field compared to its dark 
matter counterpart. The length scale below which this becomes important is the Jeans 
scale. In linear theory, the baryon overdensity obeys the following equation in a dark matter 
dominated universe (Bi et al. 1992; [Peebles 1993| ): 



+ ^H— = AnGpDMODM 4 , (Al) 



where the tilde denotes functions in Fourier space as before, H is the Hubble constant, ks 
is the Boltzmann constant, G is the Newton constant, p£,M is the average dark matter mass 
density, T is the average temperature of the gas and /i is the mean mass of each gas particle 
(for a fully ionized gas composed of hydrogen and helium with primordial abundances, it is 
about 0.6 times the proton mass). The relation between the temperature (not its average 
but its actual value) and 1 + (5b is described by 7, the temperature being proportional to 



The Jeans scale is defined in equation ([Tq) . For a dark matter dominated universe, one 
can replace pdm by the total mean density of the universe. 

For the special case of T oc a~^, making use of an equation for 5dm which is the same 
as equation ( [Al| ) except for the absence of the temperature term, it can be shown that 

r /, X 5DM(k) /.^N 

^'^^^ = JTTWW) ' ^ ^ 

if one ignores decaying modes. It expresses in a quantitative way the expectation that the 
overdensity in baryons is the same as that of dark matter on large scales (low k) but is lower 
on small scales (high k). For T with some other power-law dependence on a, solutions for 
equation ( |A1| ) are more complicated but the low and high k limits are the same: 6h = Sbm 
for small k and (5b = (5dm^j/^^ for large k (Bi et al. 1992). We note, however, that in 
practice, the gas temperature is expected to rise from almost zero before reionization to 
around 10^ — 10^ K after, and so the power-law time-dependence of T is probably not 
realized. The actual linear smoothing scale of the gas should be a little different from 
the one given above and one has to be careful about boundary conditions for 6b and its 
derivative at the onset of reionization (see Gnedin & Hui 1997). 
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B. Thermal and Ionization Evolution 



The evolution of temperature is governed by: 



-2HT 



2T d6y, 



3(1 + 5b) dt J2iX^ dt 



+ 



2 dQ 
Sksnh dt 



(Bl) 



where d/dt is the Lagrangian derivative following each fluid element, Uf, is the proper 
number density of all gas particles and T is the temperature which depends on both space 
and time. The symbol Xi is defined by = {l + 6h)Xi pb/rrip, where rii is the proper number 
density of the species i, pb is the mean mass density of baryons at the time of interest, rrip 
is the mass of the proton and (5b is the overdensity as in equation (H). The neutral fraction 
of hydrogen, Xhi (distinct from Xri) as in equation (^), is then Xm/ (-^hi + Xhii)- Note 
that Xi is a function of space and time in general. 

The first two terms on the right hand side take care of adiabatic cooling or heating. 
The third accounts for the change of internal energy due to the change in the number of 
particles. The last term dQ/dt is the heat gain (or negative heat loss) per unit volume 
by the gas particles from the surrounding radiation field. At a redshift of 2 to 4 and for 
densities of our interest, the main source of heat gain is photoionization and the main 
source of heat loss is through the recombination of ionized hydrogen and the free electron. 
At higher redshifts, other processes become important, such as Compton cooling. More 
discussion on these processes will be presented in Hui and Gnedin (1996). We note that 
one particularly simple solution of equation ( plf ) is T oc a^^(l + 5b)^^^, which holds when 
the last two terms on the right hand side can be ignored i.e. pure adiabatic expansion or 
compression. 

The above equation has to be supplemented by one that determines the abundance of 
each particle type, which takes the form: 



dXi 
IT 



—XiPh + ^ XjXkR 



:i + 6b 



nir. 



fB2l 



For instance, if Xi = Xhi, Ph is the photoionization rate. It is given by: 

Ph= 4:TcJ„am] 



dv 
'hi/ 



(B3) 



where h is the Planck constant, hz/ni = 13.6 eV, cthi is the cross-section for photoionization 
as a function of the frequency v and Ji, is the specific intensity. The photoionization rate 
Ph depends on the normalization as well as spectrum of J^^. The specific intensity Jy is 



generally taken to have a power law spectrum, v ^ io v for frequencies just above 
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z/hi- The spectrum at higher frequencies is less important for the photoionization rate of 
hydrogen. A convenient way to hide our ignorance of the spectrum is to define Jhi as in 
equation 



For Xi = Xhi, R is the recombination rate of ionized hydrogen and the free electron 
{Xj = Xf. and X^ = Xmi in equation |]B^] ): 



For Jhi with the values noted in Sec. |3.2| , the photoionization time-scale is much 
shorter than the Hubble time. This means that hydrogen is highly ionized and is essentially 
in ionization equilibrium. The two terms on the right hand side of equation (p^) almost 
balance each other, which implies equation (ITT 



We now have all the equations in place to compute the thermal and ionization 
evolution. The overdensity 5b is evolved using the Zel'dovich approximation. Its rate of 
growth is substituted into equation (pi]), which is solved together with equation (|B^ ). The 
initial conditions are as follows. The gas temperature T is equal to the cosmic microwave 
background temperature at 2; = 100 (maintained by Compton scattering) and evolves 
adiabatically after that until the universe is reionized by the UV background. Abundances 
are assumed to be primordial, which is consistent with observations so far for column 
densities less than about 10^^'^ cm~^ (See Songaila and Cowie 1996. Cooling processes due 
to metals are not important for our densities of interest in any case). All species are neutral 
until reionization occurs. One can integrate equations ( [BID and ( |B2D forward starting from 
any time between z = 100 and the beginning of reionization to obtain T and Xhi or Xhi. 
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Label 




n 


A;s/Mpc-i 








CDMla 


0.7 


1.0 


2.3 


0.05 


1.12 


-2 


CDMlb 


0.7 


1.0 


2.3 


0.06 


1.12 


-2 


CDM2 


0.8 


1.0 


1.7 


0.05 


1.13 


-1.94 


CDM3 


0.57 


1.0 


1.7 


0.05 


0.81 


-1.94 


CDM4 


0.46 


1.0 


4.9 


0.05 


1.01 


-2.17 


CDM5 


0.8 


0.5 


8.4 


0.05 


1.01 


-2.59 



Table 1: A list of all the CDM models discussed in this paper. All have h = 0.5. The larger 
scale power-spectral index is n. Every model has Qf, = 0.05 except for CDMlb, which has a 
higher baryon content and is shown in Fig. |^. The truncation scales ks for CDMla, CDMlb 
and CDM2 are defined by = 1.5 /cnl (eq. The smoothing scales for the rest of the 

models are chosen for the purpose of testing (see Sec. 0) but always larger than the length 
scale defined above. The rms density fluctuation ctq and the small scale effective slope of the 
power spectrum Ues are defined in equations (|37D and (^), for the value of ks given here. 
We assume z = 3. The rms linear density fluctuation in a sphere of radius 8h~^ Mpc is equal 

to (Tg. 
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Label 




h 


n 


Qrms/ /^K 


T/S 




0-8 




^cff 


Al 


0.1 


0.5 


0.95 


18.5 


7(1 -n) 


3.8 


0.73 


1.08 


-2.23 


A2 


0.2 


0.5 


0.95 


18.5 


7(1 - n) 


8.4 


0.66 


0.84 


-2.48 


A3 


0.3 


0.5 


0.95 


18.5 


7(1 -n) 


not apply 


0.64 


not apply 


not apply 


Bl 


0.1 


0.5 


0.9 


19.2 





3.3 


0.82 


1.08 


-2.25 


B2 


0.2 


0.5 


0.9 


19.2 





8.4 


0.73 


0.87 


-2.52 


CI 


0.1 


0.65 


0.9 


19.2 


7(1 - n) 


3.1 


0.82 


1.09 


-2.14 


C2 


0.2 


0.65 


0.9 


19.2 


7(1 -n) 


10.9 


0.76 


1.02 


-2.5 


Dl 


0.1 


0.65 


0.8 


20.5 





2.8 


0.94 


1.08 


-2.19 


D2 


0.2 


0.65 


0.8 


20.5 





10.9 


0.87 


1.02 


-2.57 



Table 2: A list of all the CHDM models discussed in this paper. All have Qb = 0.05. The 
parameters are defined as follows: Qi, is the density parameter in neutrino, n is the large 
scale spectral index of the power spectrum, Qrms is the COBE quadrupole in and T/ 5* is 
the tensor to scalar ratio. The smoothing wavenumber ks for each Q,^ = 0.1 model is 1.5 /cnl 
(eq. [|I^) and for each = 0.2 model is the Jeans scale for 7 = 1.5, Tq = 10^ and 
the corresponding h (see Sec. |3.2| ). No simulation is run for A3, so no ks is listed. The rms 
linear density fiuctuation ctq and the small scale effective spectral slope rics are defined in 
equations (|37D and (p9D , for each value of k^ given in the seventh column. The models and 
their power spectra are taken from Ma (1996). 
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Fig. 1. — A line of sight through a erg = 0.7 CDM simulation (produced using the truncated 
Zel'dovich approximation) at ^ = 3, with h = 0.5. The transfer function is taken from 
Ma (1996). Box size is 12.8 Mpc with a grid spacing of 0.032 Mpc. The parameters are 
flbh^ = 0.0125, Jhi = 0.5, Tq = 10^ K ,7 = 1.5 and fcg = 2.3 Mpc^^. All distances are 
comoving. See Sec. |^ for definitions of the symbols. The abscissas for the lower two panels 
are the comoving distances along the line of sight in units of Mpc. The lower of the two 
panels is the profile of overdensity 6h (eq. ^) and the upper one is the profile of velocity u 
(eq. 0]). The top two panels are both transmission profiles where r is the Lya optical depth 
and the abscissas represent Uo, which is related to the observed frequency through equation 
(|). The profile with solid line is obtained using the full density and peculiar velocity fields. 
The profile with dashed line is obtained using the same density field but setting the peculiar 
velocity to zero everywhere (in which case, u becomes linear in x). 
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Fig. 2. — The same as in Fig. |1| but a different line of sight. 
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Fig. 3. — The column density distributions for the same model and parameters as those in 
Fig. |1|. The quantity (P Ni^ya/ dNm/ dz has units of cm^. Crosses represent the distribution 
obtained by applying the Threshold-Algorithm (threshold at 0.89, which is the mean 
transmission) to spectra generated using the density and peculiar velocity fields predicted 
by the truncated Zel'dovich approximation. Open triangles represent the same except that 
peculiar velocities are set to zero. Open squares are obtained by applying the Threshold- 
Deblending-Algorithm at the threshold of 0.89. 
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Fig. 4. — Same parameters as in Fig. |l]. The column densities computed using two different 
methods are plotted against each other. First, we identify absorption lines by the Threshold- 
Deblending-Algorithm using a transmission threshold of 0.89 and assign column densities 
according to equation (PT|), which are plotted as the abscissas. We then take each absorption 
line identified using the Threshold-Deblending- Algorithm and search for the corresponding 



maximum in (5b and apply equation (|27|) to assign a second set of column densities, which 
are plotted as the ordinates. 
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Fig. 5.— The model CDMlb (Table ^ at z = 3, Jm = 0.325, Tq = lO^K and 7 = 1.45. 
Solid triangles represent the distribution obtained by applying the Voigt-profile-fitting- 
technique to synthetic spectra from a full hydrodynamic simulation (Zhang et al. 1996) 
with box size of 9.6 Mpc comoving and grid spacing of 0.075 Mpc. Open triangles and open 
squares are the predictions of the Density-Peak-Ansatz (DPA) coupled with the truncated 
(fcs = 2.3 Mpc~^) Zel'dovich approximation, the former using the same box size and grid- 
spacing as the hydrodynamic simulation and the latter using a box size of 12.8 Mpc and 
grid spacing of 0.05 Mpc. Crosses are the results of applying the Density-Peak-Ansatz and 
the Threshold-Algorithm (TA; transmission threshold at 0.83, the mean transmission) to the 
same density field as for the open squares. The short-dashed and long-dashed curves are 
the predictions of the Density-Peak-Ansatz coupled with the lognormal approximation, the 
former with ks = 2.3 Mpc~^ and the latter with the smoothing scale chosen so that the final 
rms density fluctuation matches that of the Zel'dovich approximation [ks = 3.6Mpc~^). 
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Fig. 6. — The solid triangles are same as those in Fig. representing the column 
density distribution obtained by applying the Voigt-profile-fitting-technique to synthetic 
spectra from a full hydrodynamic simulation. The rest of the points represent the 
column density distributions obtained using the Density-Peak-Anatz in conjunction with 
the truncated Zel'dovich approximation for two simulations of different resolutions but the 
same cosmological, thermal and ionization parameters as in Fig. ^. The open squares and 
open triangles represent the distributions of a simulation with box size of 12.8 Mpc and grid 
spacing of 0.05 Mpc (all distances quoted are comoving). The open squares here are the 
same as those in Fig. ^ where each maximum over three cells is identified as a peak. The 
open triangles are the result of a different definition of peaks: a local maximum over five 
cells with the density slope on each side of the maximum not changing signs. Similarly, the 
crosses and open hexagons are the distributions for a simulation of box size 12.8 Mpc and 
grid spacing 0.0284 Mpc, using the three-cell and five-cell definition of peaks respectively. 
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Fig. 7. — The column density distribution of the as = 0.7 CDM model with no tilt 
(CDMla in Table P, obtained using the Density-Peak- Ansatz and the truncated Zel'dovich 
approximation {ks = 2.3Mpc~^). The redshift is z = 3. Box size is 12.8 Mpc with grid 
spacing of 0.05 Mpc. Open squares: F = 0.25 (eq. [^). Crosses: F = 1. Open triangles: 
F = 5. All of them have the equation of state described by 7 = 1.5 (eq. [l^). The points 
with error bars are the observational data at about z = 3 which have been corrected for 
incompleteness, taken from Hu et al. (1995). 
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Fig. 8. — Column density distributions of the CDMla model (Table |1]) for three different 
equations of state (eq. |TP|). Open squares: 7 = 1.2; crosses (same as crosses in Fig. |^): 
7 = 1.5; open triangles: 7 = 1.7. F = 1 (defined in eq. ||32|) for all three. Points with 
error-bars are the same observational data as in Fig. |^. Long-dashed and short-dashed lines 
have the approximate slopes {f5 = 1.62 and (3 = 1.48) (normalization is chosen by hand) as 
given in equation (^) for the open triangles and open squares respectively. 
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Fig. 9. — The effect of the amphtude of power spectrum on the slope of the column density 
distribution. Open trianges represent the CDM2 model with ctq = 1.13 and ries = —1.94 
and open squares represent the CDM3 model with ctq = 0.81 and the same Ues (see Table |I|). 
The normalizations of the column density distributions are chosen by hand. The solid and 
dashed lines have the approximate slopes (/5 = 1.53 and P = 1.72) as given in equation (^Tj) 
for the open triangles and squares respectively. 
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Fig. 10. — The effect of the shape of the power spectrum on the slope of the column density 
distribution. Open trianges represent the CDM4 model with ctq = 1.01 and ries = —2.17 and 
open squares represent the CDM5 model with the same (Tq and Ues = —2.59 (see Table |1|). 
The normalizations of the column density distributions are chosen by hand. The solid and 
dashed lines have the approximate slopes (/5 = 1.67 and P = 1.79) as given in equation (^Tj) 
for the open triangles and squares respectively. 
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Fig. 11. — Column density distributions of the CDMla model (Table [1|) for three different 
initial smoothing scales. Crosses (same as crosses in Fig. ks = 2.3 Mpc~^, which is 



the smoothing scale according to the standard prescription (eq. [15|)- Open triangles: 
ks = 1.15 Mpc-^ Open squares: k^ = 8.4 Mpc ^, which is the Jeans scale for Tq = 1 and 
7 = 1.5. Significant amount of orbit-crossing has probably taken place for the last case. We 
adopt F = 1 (eq. [Q) and 7 = 1.5 (eq. [0) for all three cases. 
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Fig. 12. — (To versus ks (eq. fS^]) at z = 3. Going from the top, the sohd hne is the same 
as = 0.7, h = 0.5 CDM model with no tih as in Fig. the next four sets of points/hnes 
close together are all Q,^ = 0.1 CHDM models, Al, Bl, CI and Dl in Table |^ ; the next four 
sets are all fl^ = 0.2 CHDM models, A2, B2, C2 and D2; the last set, solid triangles, is an 
= 0.3 CHDM model, A3. Note that ao{ks) oc (gee eq. [H). 
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Fig. 13. — Column density distributions of the B2 CHDM model (see Table ^) for three initial 
smoothing scales. Crosses: ks = 19.2 Mpc~^ (standard truncation prescription, ks = 1.5 /cnl 
according to eq. ||l5l). Open triangles: k^ = 8.4Mpc~^ (Jeans scale for Tq = 10^ and 
7 = 1.5). Open squares: no smoothing at all. Points with error-bars are the observational 
data as in Fig. ^. For all models, F = 1 (eq. |^) and 7 = 1.5 (eq. [l^) are used. 
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Fig. 14. — Column density distributions for four Qi, = 0.1 CHDM models. Points with 
error-bars are the observational data as in Fig. |^. All models have Qi, = 0.1. We use 7 = 1.5 
in the equation of state for all of them (eq. |]19|)- Table |^ contains descriptions of each of the 
following models. Open hexagons: Bl, F = 3.33. Open triangles: Al, F = 3.33. Crosses: 
Dl, F = 5.7. Open squares: CI, F = 5.7. F is defined in equation (p2D. 
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Fig. 15. — The column density distribution for the A2 CHDM model (Table Three 
values of F (eq. [^]) are shown: F = 1 (open squares), F = 2.5 (crosses) and F = 5 
(open triangles). We choose 7 = 1.5 for all three (eq. Il9|). Points with error-bars are 
the observational data as in Fig. |^. The dashed line has the slope of (3 = 1.86, as given in 
equation (|4l|) . The normalization of the line is chosen by hand. 
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Fig. 16. — The column density distribution for the B2 CHDM model (Table Three 
values of F are shown (eq. |3^): F = 1 (open squares), F = 2.5 (crosses) and F = 5 (open 
triangles). 7 = 1.5 for all three (eq. [0). Points with error-bars are the observational data 
as in Fig. ^. The dashed hne has a slope of /3 = 1.86, as given in equation (^TJ). 
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Fig. 17. — The column density distribution for the D2 (Table ^ CHDM model. Three values 
of F are shown (eq. |^): F = 7.14 (open triangles), F = 3.57 (crosses) and F = 1.79 (open 
squares). 7 = 1.5 for all three (eq. [0). Points with error-bars are the same observational 
data as in Fig. 0. The dashed line has a slope of /? = 1.78, as given in equation (^1]). 
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Fig. 18. — The column density distribution of the A2 CHDM model (Table ^ for three 
different values of 7 (eq. [|l^]). F = 2.5 (eq. [^) is adopted. Three values of 7 are shown: 
7 = 1.2 (open squares), 7 = 1.5 (crosses) and 7 = 1.7 (open triangles). Points with error- 
bars are the observational data as in Fig. ^. The dashed line has a slope of (3 = 1.75, which 
is the value given in equation ( p|) for 7 = 1.2 and ctq as given in Table 



